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Abstract 

We derive Painleve-type expressions for the distribution of the m*^ largest eigen- 
value in the Gaussian Orthogonal and Symplectic Ensembles in the edge scaling limit. 
This work generalizes to general m the m = 1 results of Tracy and Widom [23] . The 
results of Johnstone and Soshnikov (see ^3], PH]) imply the immediate relevance of 
our formulas for the m*^ largest eigenvalue of the appropriate Wishart distribution. 
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Chapter 1 
Introduction 



1.1 Motivation 

The Gaussian /3-ensembles are probability spaces on A^-tuples of random variables 
{ii, . . . , ijy}, with joint density functions given by * 



p;3(^i,...,M=^r(o 



Cf^ exp 



N 



L J= 



j<k 



[1.1.1) 



The ) are normalization constants, given by 



TV 



(TV) 
13 



vr 



-TV/2 2-A^-/3TV(TV-l)/4 _ JJ' 



r(i + 7)r(i + f- 
r(i + fj) 



(1.1.2) 



*Iii many places in the Random Matrix Theory Hterature, the parameter (3 (times 1/2) appears 
in front of the summation inside the exponential factor, in addition to being the power of the 
Vandermonde determinant. That convention originated in |16j . and was justified by the alternative 
physical and very useful interpretation of Hl.l.l() as a one-dimensional Coulomb gas model. In 
that language the potential W = - HiKjMh - lj\ and P^^\t) = Cexp{-W/kT) = 

C exp(— /3 W), so that (3 ^ [k T)^^ plays the role of inverse temperature. However, by an appropriate 
choice of specialization in Selberg's integral (see Section rA.3|) . it is possible to remove the (} in the 
exponential weight, at the cost of redefining the normalization constant Cp^\ We choose the latter 
convention in this work since we will not need the Coulomb gas analogy. Moreover, with computer 
simulations and statistical applications in mind, this will in our opinion make later choices of standard 
deviations, renormalizations, and scalings a lot more transparent. It also allows us to dispose of the 
pesky square root of 2 that is often present in the definition of the Tracy-Widom /3 = 4 distribution 
in the literature. More on that later. 
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By setting f3 = 1,2,4 we recover the (finite A^) Gaussian Orthogonal Ensemble 
(GOEtv), Gaussian Unitary Ensemble (GUEat), and Gaussian Symplectic Ensem- 
ble (GSEat), respectively. We restrict ourselves to those three cases in this thesis, and 
refer the reader to [H] for recent results on the general (3 case. Originally the ij are 
eigenvalues of randomly chosen matrices from corresponding matrix ensembles (see 
Section so we will henceforth refer to them as eigenvalues. With the eigenvalues 
ordered so that ij > ij+i, define 

2-1/2 jV"V6' 1,-L.i.o; 

to be the rescaled m*^ eigenvalue measured from edge of spectrum. A standard result 
of Random Matrix Theory about the distribution of the largest eigenvalue in the 
/3-ensembles (proved only in the f] = 1,2,4 cases) is that 



whose law is given by the Tracy- Widom distributions. 
Theorem 1.1.1 (Tracy, Widom [12], 111!). 



F2{s) := P^,J, <s)= exp 



{x — s) q^{x)dx 



(1.1.5) 



Fl{s) 



F2 ■ exp 



q{x)dx 



;i.i.6) 



F2 ■ cosh 



/"OO 

'2 



q{x)dx 



[1.1.7) 
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The function q is the unique (see solution to the Painleve II equation 



q = xq + 2q 



Xl.8) 



such that q{x) ~ Ai(a;) as x ^ oo, where Ai{x) is the solution to the Airy equation 
which decays like ivr^^/^x^^/^ exp (— at +oo. The density functions corre- 
sponding to the F/3 are graphed in Figure ITTTl ^ Let F2{s^m) denote the distribution 




Figure 1.1: Tracy- Widom Density Functions 

for the m*'* largest eigenvalue in GUE. Tracy and Widom showed in [22] that if we 
define -F2(s, 0) = 0, then 



Fofs, m + 1) - Fois.m) 



m! c/A™ 



D2{s,X)\^^^, m > 0, 



(1.1.9) 



where 



-^2(5, A) = exp 



{x — s) q^{x, X)dx 



(1.1.10) 



and q{x, A) is the solution to (jl.l.Sj) such that g(x, A) ~ \fX Ai(x) as x ^ 00. An 
intermediate step leading to p.l.lO|) is to first show that (s, A) can be expressed 



^The square root of 2 in the argument of reflects a normalization chosen in Ijl. 1.1(1 to agree 
with Mehta's original one (see See first footnote in this Introduction, as well as comments 

following Equation l|4.3.41|) . and Section IX4)l . 
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D2{s,X) = det{I - XKm), (1.1.11) 

where K^i is the integral operator on oo) with kernel 

^„(,,,),AiMAi-fo)-Ai-MAifa)^ (1,1,12) 

x-y 

In the /3 = 1,4 result similar to (|l.l.lip holds with the difference that the 

operators in D^(s,A) have matrix-valued kernels (see e.g. [IE])- In fact, the same 
combinatorial argument used to obtain the recurrence (jl.l.9|) in the f3 = 2 case also 
works for the (3 = 1,4 cases, leading to 

Fi3{s,m + l)-Fp{s,m)=^-^j^D'J\s,X)\^^^, m>0,/3 = l,4, (1.1.13) 

where F^(s,0) = 0. Given the similarity in the arguments up to this point and 
comparing ()1.1.10|) to ()1.1.5|) . it is natural to conjecture that -D^(s, A),/5 = 1,4, can 
be obtained simply by replacing q{x) by q{x,X) in (jl.l.fij) and ()1.1.7j) . However the 
following fact, of which Corollary ()1.2.4|1 gives a new operator theoretic proof, hints 
that this cannot be the case. 

Theorem 1.1.2 (Baik, Rains j2]). the appropriate scaling limit, the distribution 
of the largest eigenvalue in GSE corresponds to that of the second largest in GOE. 
More generally, the joint distribution of every second eigenvalue in the GOE coincides 
with the joint distribution of all the eigenvalues in the GSE, with an appropriate 
number of eigenvalues. 

This so-called "interlacing property" between GOE and GSE had long been in the 
literature, and had in fact been noticed by Mehta and Dyson (see ^Zj). In this context, 
the remarkable work of Forrester and Rains in ^2] classified all weight functions for 
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which alternate eigenvalues taken from an Orthogonal Ensemble form a corresponding 
Symplectic Ensemble, and similarly those for which alternate eigenvalues taken from 
a union of two Orthogonal Ensembles form a Unitary Ensemble. Theorem ()1.1.2|1 does 
not agree with the formulae we postulated for Di3{s,X),j3 = 1,4. Indeed, combining 
the two leads to incorrect relationships between partial derivatives of q{x, A) with 
respect to A evaluated at A = 1 and q{x, 1). To be precise, the conjecture is true for 
D^^s, A) but it is false for -Di(s, A). The correct forms for both Dp{s, X),(3 = 1,4 are 
given below in Theorem ()1.2.H) . 

This work also extends that of Johnstone in ^S] (see also JO]), since Fi{s,m) 
gives the asymptotic behavior of the m*'' largest eigenvalue of a p variate Wishart 
distribution on n degrees of freedom with identity covariance. This holds under very 
general conditions on the underlying distribution of matrix entries by Soshnikov's 
universality theorem (see ^Hj for a precise statement). In Table EiSl we compare our 
distributions to finite n and p empirical Wishart distributions as in [12] . 

1.2 Statement of the Main Results 

Theorem 1.2.1. In the edge scaling limit, the distributions for the tv}^ largest eigen- 
values in the GOE and GSE satisfy the recurrence p.l.l3p with^ 



•I- See first footnote in this Introduction, as well as comments following Equation H4.3.41|l . and 
Section I A. 41 for an explanation of why we write Di{s, A) instead of Di{s/\/2, A) as is customary in 
the RMT literature. 



Di(s,A) =D2(S,A) 



A — 1 — cosh /i(s. A) + \/a sinh /i(s. A) 



(1.2.1) 



A-2 




(1.2.2) 
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Figure 1.2: 10^ realizations of 10^ x 10^ GOE matrices; the solid curves are, from right 
to left, the theoretical limiting densities for the first through fourth largest eigenvalue. 

where 

/oo 
q{x,\)dx, A:=2A-A^ (1.2.3) 

and q{x, A) is the solution to (|1.1.8|) such that q{x, A) ~ Ai(x) as x — ^ oo. 

Corollary 1.2.2 (Interlacing property). 

F4{s,m) = Fi{s,2m), m > 1. (1.2.4) 

In the next section we give the proof of these theorems. In the last, we present an 
efficient numerical scheme to compute F^(s, m) and the associated density functions 
fpis, m). We implemented this scheme using MATLAB™ ,^ and compared the results 
to simulated Wishart distributions. 



§MATLAB™ is a registered trademark of The MathWorks, Inc., 3 Apple Hill Drive, Natick, 
MA 01760-2098; Phone: 508-647-7000; Fax: 508-647-7001. 



Chapter 2 



Preliminaries 



2.1 Determinant matters 

We gather in this short chapter more or less classical results for further reference. 
Since these are well-known results, their proofs will be omitted. The exception to 
this is the first one which is a simple identity relating a determinant to the fourth 
power of a Vandermonde, whose proof is given because it is a little less common. 

Theorem 2.1.1. 



Yl {xj - Xk^ = det {xl j xi 



)i=0,...,2Ar-l 
k=l,...,N 



0<j<k<N 



Proof. First note that 



det {xl j xl ^) 



d 



det {x>^ yi)j=o,...,2N-i 



j=0,...,2N-l — 
k=l,...,N 



dyi... 



k=l,...,N 



yk=^k 
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Now 



det {xi yi)j=o,...,2N^i 

k=l,...,N 



1 1 1 1 ■■■ 

xi yi X2 y2 ■■■ 



A yl xl yl ■■■ 



\ 



\ 



J 



JJ [xj - Xfc) {yj - Vk) ivj - Xk) {xj - yk) JJ {yj - 

j>k 



Xj) 



Differentiate this product with respect to each yi and applying the product rule, we 
obtain 



d 



dyi . . . dyN 



W^Vj - Vk) n ~ ~ ~ y''^ n ^y^ 



j>k 



+ 



+ 



+ 



d 



dyi 


■ ■dyN 




d 


dyi 


■ ■dyN 




d 


dyi 


■ ■dyN 



Yliyj - Xk) n (^j' ~ ^'^^ ~ y^^ ~ y''^ W ^y^ ~ 



\\{xj -yk) W {xj - Xk) {yj - yk) {y^ ~ Xk) iVj - xj) 



W^Vi - ^k) n ~ ^'^^ ^yj ~ y^^ ^y^ " ^^'> ~ y^^- 

j>k 



We set yk = Xk after differentiating. In each of the first three terms the last factor 
wiU be 0. Therefore, only the last term will remain with yk = Xfc. Now 



d 



dyi... dyN 



Yliyj - Xk) 



1. 



So we are left with 



JJ {xj - Xk) (yj - yk) iVj - Xk) (xj - yk) 

j>k 



yk=^k 



Y[ (Xj - Xk)^. 

0<j<k<N 



□ 
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Theorem 2.1.2. If A,B are Hilhert-Schmidt operators on a general * Hilhert space 
H, then 

det{I + AB) = det{I + BA). 
Theorem 2.1.3 (de Bruijn, 1955). 



det((^j(xfe))^^^ det(V'j(a;fc)),<^, c?/i(a;i) ■ ■ ■ c?/i(x7v) 
= A^! det I / (pj{x) ipk{x) d fi{x) 



(2.1.1) 



l<j,k<N 



Xl<...<Xj^ 



det((^j(xfc))^<^ ,<^d/i(xi) ■ ■ -dfi^XN) 

= Pf( / / sgn{x -y)ipj{x)(pk{x)dn{x)dfi{y) 



(2.1.2) 



l<j,k<N 



det{(pj{xk) ipj{xk))^^.^2j^djj,{xi) ■ ■■d^{xN) 

r<k<N 



(2 AT)! Pf / (^,(x)V^fc(x) -<^fc(x)V,(a:)rf/i(a:; 



(2.1.3) 



l<j,fc<2 iV 



where Pf denotes the Pfaffian. The last two integral identities were discovered by 

de Bruijn jB] in an attempt to generalize the first one. The first and last are valid 

in general measure spaces. In the second identity, the space needs to be ordered. In 

the last identity, the left hand side determinant is a 2 x 2 determinant whose 

columns are alternating columns of the (pj and ifjj (i.e. the first four columns are 

{ifj{xi)], {^/'j(xi)}, {(pj{x2)], {%l)j{x2)], respectively for j = 1,...,2A^), hence the 

notation, and asymmetry in indexing. 
*See [ni for proof. 
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2.2 Recursion formula for the eigenvalue distribu- 
tions 



With the joint density function defined as in (jl.l.lj) . let J denote the interval (t, oo), 
and X = Xj{x) its characteristic function. "I" We denote hj X = 1—X the characteristic 
function of the complement of J, and define X^ = 1 — XX. Furthermore, let E/3^N(t, m) 
equal the probability that exactly the m largest eigenvalues of a matrix chosen at 
random from a (finite N) /3-ensemble lie in J. We also define 

Gf3,N{t,\) = j " J X^{Xi) ■ ■ -X^ixN) Pf3{Xi, . . . ,X]y) dxi ■ ■ ■ dxN- (2.2.1) 

For A = 1 this is just -E/3,7v(t, 0), the probability that no eigenvalues lie in (t, oo), or 
equivalently the probability that the largest eigenvalue is less than t. In fact we will 
see in the following propositions that Gp^iqit, A) is in some sense a generating function 
for Ep^^it, ni). 

Proposition 2.2.1. 

G'ftiv(t,A) = J](-A)'=(^^^ j ■■■ jp^{xu...,XM)dxi---dxN. (2.2.2) 

Proof. Using the definition of the X^(xi) and multiplying out the integrand of 1)2.2.11) 
gives 

Gf3,N{t, A) = ^(-A)^'y" ■■■ j ek{X{xi), X{xn))P(3{xi, . . . , xn) dxi---dxN, 
where, in the notation of Ck = rriik is the /c*'* elementary symmetric function. 



^Much of what is said here is stiU vahd if J is taken to be a finite union of open intervals in 
M (see 1^). However, since we wiU only be interested in edge eigenvalues we restrict ourselves to 
(t, oo) from here on. 
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Indeed each term in the summation arises from picking k of the A X-terms, each of 
which comes with a negative sign, and N — k of the I's. This explains the coeffi- 
cient (— A)'^. Moreover, it follows that contains (^) terms. Now the integrand is 
symmetric under permutations of the Xi. Also if Xj ^ J, all corresponding terms in 
the symmetric function are 0, and they are 1 otherwise. Therefore we can restrict 
the integration to Xj G J, remove the characteristic functions (hence the symmetric 
function), and introduce the binomial coefficient to account for the identical terms 
up to permutation. □ 

Proposition 2.2.2. 



Ep^N{t,m) 



-1' 



ml dX^ 



Gl3,N(t, A) 



m > 0. 



(2.2.3) 



A=l 



Proof. This is proved by induction. As noted above, Ej3^N(t,0) = G/3,Ar(t, 1) so it 
holds for the degenerate case m = 0. When m = 1 we have 



— ^ Gp^N{t, A) 



A=l 



d 

"dX 

N 

E- 



X^(.Xi) ■ ■ -X^iXn) Pp^\x) dXi - --dXri 

X(xi)---X(x,_i)X(x,)X(x,+i)--- 



A=l 



■■■X{xn)Pp (x) dxi ■ ■ ■ dxj\f. 

The integrand is symmetric under permutations so we can make all terms look the 
same. There are N = (^) of them so we get 



-^G,,^.(t,A) 



A=l 



A^' 



X{xi) X{x2) ■ ■ ■ X{xn) Plf^\x) dxi - ■ ■ dxN 
X{xi) X{x2) ■ ■ ■ X{xn) Pp^\x) dxi - ■ ■ dxN 



A=l 



E/3,N{t, 1) 
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When m = 2 then 



dX 



N 



N 

i=2 



■ ■ ■ X{xn) p'i^^\x) dxi - ■ ■ dxN 
N{N -1) 



A=l 



N 
2 

N 



X{xi) X{x2) X{x3) ■ ■ ■ X{xn) Pp^\x) dxi - ■ ■ dxM 
X{xi) X{x2) X{x^) ■ ■ ■ X{xn) P^fl^\x) dxi - ■ ■ dxN 
X{xi) X{x2) X{x3) ■ ■ ■ X{xn) P^^\x) dxi - ■ ■ dxN 



A=l 



A=l 



Ei3,N(t, 2), 



where we used the previous case to get the first equahty, and again the invariance of 
the integrand under symmetry to get the second equahty. By induction then, 



1 

m ! 



d 



Gp,Nit, A) I 



A=l 



N{N -1)---{N -m. + 2) 



m I 



N 

E 

j=ni 



X{xi)X{x2) 



■ ■ ■ X{xj_i) X{xj) X{xj+i) ■ ■ ■ X{xn) Pf3^\x) dxi - ■ ■ dxN 
N{N -1)---{N -m+1) 



m I 



X{xi) 



■X{Xm) X{Xm+l) ■ ■ ■X(xAr) P'p^\x) d Xi ■■■dXN 



N 



m 



X{xi) ■ ■ ■ ■ ■ ■ X{xn) Pf\x) dxi---dxN 



Ef3,N(,t,m). 



□ 

If we define Fji^]^{t,m) to be the distribution of the im}^ largest eigenvalue in the 
(finite A^) /^-ensemble, then the foHowing probabihstic resuh is immediate from our 
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definition of Ef^^j^(t,m). 
Corollary 2.2.3. 

Fp^^it, m + 1) - m) = E^^^it, m) (2.2.4) 



Chapter 3 



The distribution of the m largest 
eigenvalue in the GUE 

3.1 The distribution function as a Fredholm deter- 
minant 



We closely follow for the derivations which follow. The GUE case corresponds to 
the specialization /3 = 2 in (ll.l.lj) so that 



where w{x) = exp (— x^), f{x) = —X Xj{x), and C2 depends only on A^. In the steps 
that follow, additional constants depending solely on (such as A^!) which appear will 
be lumped into C2^\ A probability argument will show that the resulting constant 
at the end of all calculations simply equals 1. Expressing the Vandermonde as a 
determinant 



G2,N{'t, A) — 6*2 




(3.1.1) 



n -Xk) = det {xi 




(3.1.2) 



^<j<k<N 



14 
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and using ()2.1.1|) with (pj{x) = ipj{x) = and dfi{x) = w{x){l + /(x)) yields 

G2,7v(t,A) = ci^^det ( I x^^^w{x){l + f{x)) dx)] . (3.1.3) 

Let {^j{x)"\ be the sequence obtained by orthonormahzing the sequence {a;-' 
It follows that 

G2,7v(t, A) = Cf ) det ( I ip,{x) ipk{x) (1 + /(x)) dx)\ (3.1.4) 

= Cf)det(5,-fc+ /<^,(x)<^fc(x)/(x)dx)') . (3.1.5) 

The last expression is of the form det (J + AB) for A : L^(]R) with kernel 

A{j,x) = ipj{x)f{x) whereas B : — L^(M) with kernel B{x,j) = ^Pj{x). Note 
that AB -.C^ ^ has kernel 

AB{j,k)= ! ipj{x)M^)fix)dx (3.1.6) 
Jr 

whereas BA : L'^(R) L'^{M) has kernel 

N-l 

BA{x, y) = ^ V5fc(x) (pk{y) ■■= K2,n{x, y). (3.1.7) 

From lTOl it follows that 

G'2,7v(t, A) = Cf) det (/ - K2,N f) , (3.1.8) 

where K2^n has kernel K2,N{x,y) and K2^n f acts on a function by first multiplying 
it by / and acting on the product with K2^n- From ()3.1.1|) we see that setting / = 
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in the last identity yields Cg^'' = 1. Thus the above simplifies to 



16 



G2,N{t,^) = det{I-K2,Nf). (3.1.9) 

3.2 Edge scaling and differential equations 

Following the derivation in [2^1, we specialize w{x) = exp(— x^), f{x) = —XXj{x), 
so that the {Lpj{x)} are in fact the Hermite polynomials times the square root of the 
weight. Using the asymptotics of Hermite polynomials, it follows that in the so-called 

edge scaling limit, 



(3.2.1) 

is KA_i{x,y) as defined in (jl.l.l2|l . The fact that this limit holds in trace class norm 
is crucial (see jS], for proofs, and [22] ioT corresponding results in the /5 = 1,4 
cases). It is what allows us to take the limit inside the determinant in ()3.1.9p . For 
notational convenience, we denote the corresponding operator K\i by K in the rest 
of this subsection. We also think of K as the integral operator with kernel 

x-y 

where ip{x) = a/A Ai(x), ip{x) = Ai'(x) and J is (s, oo) with 

t = V2N + , . (3.2.3) 

Note that although K{x, y), ip and ip are functions of A as well, this dependence will 
not affect our calculations in what follows. Thus we omit it to avoid cumbersome 



CHAPTER 3. THE M^" LARGEST EIGENVALUE IN THE GUE 17 
notation. The Airy equation implies that (f and ip satisfy the relations 

d 

a X 

= xif. (3.2.4) 

ax 

We define D2^n{s, A) to be the Fredholm determinant det (J — K). Thus in the edge 
scaling limit 

E2{S,X) :=E2,oo(s,A) = D2(S,A). 

We define the operator 

R={I-Ky^K, (3.2.5) 

whose kernel we denote R{x, y). Incidentally, we shall use the notation = in reference 
to an operator to mean "has kernel". For example R = R{x, y). We also let M stand 
for the operator whose action is multiplication by x. It is well known that 



^logdet(J-i^) = -/2(s,s). (3.2.6) 
as 



For functions / and g, we write f ® g to denote the operator specified by 

f^g = fix)g{y), (3.2.7) 

and define 



Qix,s) = Qix) = {{I-Kr'^)ix), 
P{x,s) = P{x)= {{I - K)-^%l)){x). 



(3.2.8) 
(3.2.9) 
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Then straightforward computation yields the following facts 

[M,K] = ip tp - tp ® ifi, 
[M ,{I - Ky^] = {I - Ky^[M ,K]{I - K)~^ 

= Q®P-P®Q. (3.2.10) 

On the other hand if (/ — K)~^ = p{x, y), then 

p{x,y) = 5{x-y) + Rix,y), (3.2.11) 

and it follows that 

[M,{I-Ky']={x- y)p{x, y) = {x- y)R{x, y). (3.2.12) 

Equating the two representation for the kernel of [ M , (/ — Ky^ ] yields 

, , Q{x)P{y) - P(x)Qiy) 
R(x,y) = ' wvvy;_ 3.2.13 

x-y 

Taking the limit y — > x and defining q{s) = Q{s, s), p{s) = P{s, s), we obtain 

R{s, s) = Q'{s, s) p{s) - P'{s, s) q{s). (3.2.14) 
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Let us now derive expressions for Q'{x) and P'{x). If we let the operator D stand for 
differentiation with respect to x, 

Q'{x,s) = D{I-Ky\ 

= {I -Ky'D^+[D,{I -K)-^]^ 

= {i-Kr'i, + [D,{i-Kr']^ 

= P{x)+[D,{I-Ky^]^. (3.2.15) 

We need the commutator 

[D, (/ - KT'] ={I- KY' [D, K] (/ - KY' . (3.2.16) 
Integration by parts shows 

[D, K] = {^ + ^^+ K{x, s)5{y - s). (3.2.17) 
The 6 function comes from differentiating the characteristic function x- Moreover, 

Thus 

[D, (/ - KY'] = -Q{x)Q{y) + R{x, s)p{s, y). (3.2.19) 
(Recall (/ — KY^ = p{x, y).) We now use this in ()3.2.15j) to obtain 



Q'{x,s) = P{x)-Q{x){Q,^) + R{x,s)q{s) 
= P{x) - Q{xY{s) + R{x,s)q{s), 
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where the inner product {Q, ip) is denoted by u{s). Evaluating at x = s gives 
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Q'{s, s) = p{s) - q{s)u{s) + R{s, s)q{s). 



(3.2.20) 



We now apply the same procedure to compute P'. 



P'{x,s) = D{I-Ky^i) 

= {I-K)-^Di>+[D,{I-K)-^]^ 

= M{I- K)~' ^+[{1 - Ky\M]^+ [D, (/ - K)-'] ^ 

= xQ{x) + {P®Q-Q®P)^ + {-Q ® Q)tp + R{x, s)p{s) 

= xQ{x) + P{x) (Q, ^) - Q{x) (P, ^) - Q{x) (Q, 4j) + R{x, s)p{s) 

= xQ{x) - 2Q{x)v{s) + P{x)u{s) + R{x, s)p{s). 



Here v 



(P, 9?) = {ip, Q). Setting a; = s we obtain 



P'(s, s) = sq{s) + 2q{s)v{s) + p{s)u{s) + R{s, s)p{s). 



(3.2.21) 



Using this and the expression for Q'{ 



s, s) in ()3.2.14p gives 



R{s, s) = p^ — sq^ + 2(fv — 2pqu. 



(3.2.22) 



Using the chain rule, we have 




(3.2.23) 
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The first term is known. The partial with respect to s is 



dQ{x, s] 



{I - KY^ ^ [I - Ky\ 



ds ds 
= -R{x,s)q{s) 



where we used the fact that 

dK 



-K{x,s)S{y - s). (3.2.24) 



ds 

Adding the two partial derivatives and evaluating a.t x = s gives 



A similar calculation gives 



^S.=p-qu. (3.2.25) 

a s 



^ = sq-2qv+pu. (3.2.26) 
a s 



We derive first order differential equations for u and v by differentiating the inner 
products. Recall that 

(p{x)Q{x, s) dx. 

Thus 



— = -^{s}q{s) + J ip{x) — dx 



- \^ip{s) + J R{s,x)(p{x) dx ] q{s) 

- (I - K)-' ^(s) qis) 



Similarly, 



dv 

— = -pq. 
d s 



(3.2.27) 



CHAPTER 3. THE M^" LARGEST EIGENVALUE IN THE GUE 22 

From the first order differential equations for q, u and v it follows immediately that 
the derivative of — 2v — is zero. Examining the behavior near s = oo to check 
that the constant of integration is zero then gives 

-2v = q^. (3.2.28) 

We now differentiate ()H.2.25j) with respect to s, use the first order differential equations 
for p and m, and then the first integral to deduce that q satisfies the Painleve II 
equation 

q" = sq + 2gl (3.2.29) 

Checking the asymptotics of the Fredholm determinant det(/ — K) for large s shows 
we want the solution to the Painleve II equation with boundary condition 

g(s, A) ~ v^Ai(s) as s ^ oo. (3.2.30) 



That a solution q exists and is unique follows from the representation of the Fredholm 
determinant in terms of it. Independent proofs of this, as well as the asymptotics as 
s ^ oo were given by [H] , , [Zl • Since [D ,{I — K)~^] = {d/dx + d/ dy) R{x, y), 
dTTTm says 

^ ^) " -^(^)^(^) + ^(^' ')P^'^ (3.2.31) 

In computing dQ{x, s)/ds we showed that 

I- (/ - K)-^ = -^R{x, y) = -R{x, s)p(s, y). (3.2.32) 
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Adding these two expressions, 

^ + ^ + ^ j R{x, y) = -Q{x) Q{y\ (3.2.33) 
and then evaluating aX x = y = s gives 

-^R{s,s) = -q\ (3.2.34) 

a s 

Integration (and recalhng (j3.2.6|) ) gives, 

d 1"^ 

— \ogdet{I-K) = - q\x,X)dx; (3.2.35) 
ds J, 

and hence, 



OO / ^OO \ POO 



\ogdet{I - K) = - J q\x,\)dxjdy = -j (x -s)g^(x. A) dx. (3.2.36) 

To summarize, we have shown that D2{s, A) has the Painleve representation 

D2{s,\) = exp {x - s)q^{x,\)dx^ (3.2.37) 

where q satisfies the Painleve II equation ()3.2.29|) subject to the boundary condition 



Chapter 4 



The distribution of the m largest 
eigenvalue in the GSE 



4.1 The distribution function as a Fredholm deter- 
minant 

The GSE corresponds case corresponds to the speciahzation = 4 in p.l.lj) so that 



where w{x) = exp (— x^), f{x) = —XXj{x), and C[ depends only on A^. As in the 
GUE case, we will lump into C[^^ any constants depending only on that appear 
in the derivation. A simple argument at the end will show that the final constant is 
1. These calculations more or less faithfully follow and expand on |21|. By 
Gi,N{i^ A) is given by the integral 



G4^Ni't, A) — C4 




(4.1.1) 



/ 



/ 



N N 




det{xi jxl 



)j=0,...,2N-l 
k=l,...,N 



Y\_w{^i) + f{xi))dxi ■ --(ixN 
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which, if we define ^Pj{x) = x^~^ w{x) (1 + f{x)) and ipj{x) = {j — 1) x-'"^ and use the 
hnearity of the determinant, becomes 



Gi,N{tA) = I ■■■ / det((/)j(xfc) %l){xk))i<j<2N dxi- ■ -dxN- 

J l<k<N 



Now using (j2.1.3|) . we obtain 



G'4,Af(t, A) = cf^'Vil I (pj{x)ipkX - (fkix)ipj{x) d 



:x 



l<],k<2 N 



Cf^^ Pf(^j {k-j)x^+''-^w{x){l + f{x))dx^ 
Ci^^ Pf ( /" {k-j)x^+''-^w{x){l + f{x))dx 



0<j,fc<2 JV-1 



where we let k k + 1 and j ^ j + 1 in the last line. Remembering that the square 
of a Pfaffian is a determinant, we obtain 



G%{t, A) = Cf ^ det (^jik- j) x'^"-' w{x) (1 



f{x)) dx 



0<j.k<2 JV-1 



Row operations on the matrix do not change the determinant, so we can replace {x^} 
by an arbitrary sequence {pj{x)} of polynomials of degree j obtained by adding rows 
to each other. Note that the general (j, k) element in the matrix can be written as 



X ] X-' — \ —— X-' I X'" 



dx I \dx 



Thus when we add rows to each other the polynomials we obtain will have the same 
general form (the derivatives factor). Therefore we can assume without loss of gen- 
erality that G\^{t, A) equals 



C^^^ det ( / [pj{x)p'^{x) - p'j{x)pk{x)\ w{x) (1 + f{x)) dx 



0<j,fc<2 JV-1 
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where the sequence {pj{x)} of polynomials of degree j is arbitrary. Let ipj = pjW^^'^ 
so that pj = ijjjW"^^'^. Substituting this into the above formula and simplifying, we 
obtain 



Glj,{t,X) = Cf^ det 



0<j,k<2 N~l 



Cf^ det [M + L] = Cl"^^ det[M] det [J + M"^ ■ L], 



where M, L are matrices given by 



M 



i^jix) ^'ki^) - i^k{x) iJj{x)) dx 



0<j,fc<2 JV-1 



L = i I (ipjix) ^lj'f,{x) - tjjkix) ^pjix)) f{x)dx 

0<j,k<2 N~l 

Note that det[M] is a constant which depends only on so we can absorb it into 
C[^\ Also if we denote 



27V-1 



it follows that 



0<j,k<2 JV-1 



k=0 



M"^ -N 
Let A : L^fM) x L'^i 



i)'k{x) - r]j{x) tpk{x)) f{x) dx 



0<j,k<2 JV-1 



C be the operator defined by the 2N x 2 matrix 



A{x) 



^ rio{x) -Vo{x) ^ 



?7i(xj -r]i[x) 
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Thus if 



9o{x) 



we have 



( 



e L^m X L^ 



Ag = A{x)g 



I {r]o9o - r]ogi) dx 
I {Vi9o - V'i9i) dx 



\ 



J 



G C 



2N 



Similarly we define B : C^" — > L^(M) x L^(]R) given by the 2 x 2n matrix 



Bix) = f 



iJo{x) v^i(x) 



... J 



Explicitly if 



then 



a 



Ba = B{x) ■ a 



( \ 



e C 



2N 



I 2N-1 \ 



i=0 
2N-1 



eL^ xL^ 



1=0 
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Observe that M"^ ■ L = AB : 



''2n , (n2n 



2N-1 



ABa 



E 

2N-1 r 

E 

j=0 



C^". Indeed 



\ 



J 



«0 



(M-i ■ L) 



a. 



Therefore, by dTT^ 



Gj^it, A) = Cl^^ det(/ + ■ L) = Cj"^^ det(/ + AS) = C^'^^ det(/ + SA) 



(N) 



where BA : 



(M) . From our definition of A and B it follows that 



V 



/ 



/ (Vi9o-vl9i) dx 
f^^'il {m%-ii9i) dx 

i—n V"' 



i=0 
2N-1 



j=0 

/ „2Af-l „2N-1 \ 



v 



i=0 



j=0 
2Ar-l 



j=0 



where K4^n is the integral operator with matrix kernel 



I 2N-1 2N-1 \ 



1=0 

2N-1 



i=0 
2N~1 



i=0 



i=0 



I 
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2N-1 



Recall that rij{x) = ^-jki^ki^) so that 



k=0 



Ki,N{x,y) 



I 1N-\ 2N-1 \ 



j,k=0 
2N-1 



j,k=0 
2N-1 



\ j,k=0 j,k=0 



J 



Define e to be the following integral operator 



e = e{x- y) 



if X > y, 
if X < y. 



(4.1.2) 



As before, let D denote the operator that acts by differentiation with respect to x. 
The fundamental theorem of calculus implies that D e = eD = I. We also define 



2N-1 



SN{x,y) = Y ^ji^)f^jk^k{y)- 



j,k=0 



Since M is antisymmetric, 



2N-1 



2N-1 



2N-1 



SN{y,x) = Y '^jiy)f^jk'^k{x) = -Y '^jiy)f^kj'ipk{x) = ~Y '^jiy)f^kj1p'k{x), 



j,k=0 



j,k=0 



j,k=0 



after re-indexing. Note that 



2N-1 2N-1 

eSNix,y) = Y (^D'ipj{x)iJ,jkipkiy) = Y '^ji^)l^jkipkiy), 

j,k=0 j,k=0 



and 



d 



-^SN{x,y) = Y ^ji^)f^jki^kiy)- 
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Thus we can now write succinctly 



^ SN{x,y) -j^SN{x,y) ^ 



V 



eSN{x,y) SN{y,x) 



(4.1.3) 



To summarize, we have shown that ^^^(t, A) = cl^"* det(/ — K^j^f). Setting / = 
on both sides (where the original definition of G^^Nit, A) as an integral is used on the 
left) shows that C[^^ = 1. Thus 



G4,N{'t, A) — A / D4^N{t, A) 



(4.1.4) 



where we define 



D4,7v(t,A) = det(/-K4,7v/), 



(4.1.5) 



and K^^N is the integral operator with matrix kernel ()4.1.3p . 



4.2 Gaussian specialization 



We would like to specialize the above results to the case of a Gaussian weight function 



w{x) = exp (x^) 



(4.2.1) 



and indicator function 



We want the matrix 



/(x) = -AX, , J=(t,cx)). 



M 



{'^j{x)'ip'k{x) - ipk{x) ipjix)) dx 



0<i,fc<2 N-l 
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to be the direct sum of copies of 



^ 1^ 



so that the formulas are the simplest possible, since then fij^ can only be or ±1. In 
that case M would be skew-symmetric so that = —M. In terms of the integrals 
defining the entries of M this means that we would like to have 



ip2j+l{x) —i'2kix) 



ij2k{x) -^i)2j+l{x) ) dx 



-5 



'j-.k 



and otherwise 



^i(^) -j-^kix) - ^pj{x) -^i^kix)] dx = 0. 
dx dx ' 



It is easier to treat this last case if we replace it with three non-exclusive conditions 



i'2j{.x) i)2k{,x) - ij2k{x) ^ilj2j{x) ] dx = 0, 

dx dx 



i'2j+l{x) ijj2k+l{x) - ^2k+\{x) ^ i}2j+\[x) ) = 0, 



(so when the parity is the same for j, /c, which takes care of diagonal entries, among 
others) and 

dx dx 



whenever |j — A;| > 1, which targets entries outside of the tridiagonal. Define 



<^fc(x) = -Hk{x) exp(-a;V2) for = J^^hX^ (4.2.2) 
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where the Hj. are the usual Hermite polynomials defined by the orthogonality condi- 
tion 



Then it follows that 



Now let 



/ Hj{x) Hk{x) e dx = c] 6j,k- 



iPj(x) (pk{x) dx = Sj^k- 



1 1 

This definition satisfies our earlier requirement that ipj = pj w^^'^ with w defined in 
fj4.2.H) . In particular we have in this case 



P2j+l{X) 



H2j+i {x) 



Let e as in ()4.1.2|1 . and D denote the operator that acts by differentiation with respect 
to X as before, so that D e = e D = I. It follows that 



d 



d 



^2. (x) ^2.+i (^) - i^) % (a^) 



dx 



dx 



XI 



dx 



We integrate the first term by parts and use the fact that 



J^<^^, (x) = if^ (x) 
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and also that (p . vanishes at the boundary (i.e. if . (±00) = 0) to obtain 



^2. (x) (x) - (x) (x) 



dx 



dx 



as desired. Similarly 



d d 

^2, + l i^2k (^) - ^2. i^) + l i^) 



dx 



1 
2 

-6. 



d 



d 



-^2, + l (X) ^ ^<^2.+l (^) +^^2, + , <^2. + l 



dx 



Moreover, 



P2j+l{x) 



H2j+i {x) 



is certainly an odd function, being the multiple of and odd Hermite polynomial. On 
the other hand, one easily checks that e maps odd functions to even functions on 
L2(M). Therefore 

P2j{x) = ^-j=eH2j+i (x) 
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is an even function, and it follows that 



^2fe (^) i^2j - (X) ^jj,, (X) 



dx 



P2j (x) ^ P2k (x) - P2k (x) ^ P2j (x) 



w{x) dx 



0, 



since both terms in the integrand are odd functions, and the weight function is even. 
Similarly, 



^2. + l ^2, + l (X) - + l (X) ^P,,^, (X) 



dx 



d 



d 



P2j+1 [X) — P2k+1 {X} - P2k+1 [X) — P2j+1 [X) 



w{x) dx 



0. 



Finally it is easy to see that if |j — A;| > 1 then 



dx dx 



dx = 0. 



Indeed both differentiation and the action of e can only "shift" the indices by 1. Thus 
by orthogonality of the ipj, this integral will always be 0. Hence by choosing 



1 1 

^2,+i (x) = -j= (x) , {x) = --j=e cp^^^, (x) , 



we force the matrix 



M 



0<j,k<2n-l 
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to be the direct sum of copies of 



^ 1^ 



Hence M^^ = —M where = {fij^k}j,k=o.2N-i- Moreover, with our above choice, 
fij^k = if j, k have the same parity or \j — k\ > 1, and fi2k,2j+i = ^jk = — /^2j+i,2fc for 
j, k = 0. . . . , N — 1. Therefore 



2N-1 



SN{x,y) 



'^ji^)f^jk'ipk{y) 

j,k=0 
N-l 



N-l 



j=0 



3=0 

j=Q j=0 



Recall that the Hj satisfy the differentiation formulas (see for example P], p. 280) 



H'^{x) = 2xH^{x)-H,^^{x) j = l,2. 



(4.2.3) 



H'Ax)=2jH,.,ix) j = l,2. 



Combining 14.2.21 and 14.2.31 yields 



<^'j{x) = x<^j{x) - ^ Vj+i{x). 



Similarly, from 14 . 2?^ and l¥. 2 . 41 we have 



(4.2.4) 



(4.2.5) 



(p'j{x) = -x<^j{x) + 2j ^ <^j-i{x) 



(4.2.6) 
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Combining ()4.2.5|) and ()4.2.6p . we obtain 



Let if 



and 



•^2 



. Then we can rewrite fj4.2.7|) as 



^' = A^ 



(4.2.7) 



where A = {aj^k} is the infinite antisymmetric tridiagonal matrix with ajj-i 



Hence, 



fc>0 



Moreover, using the fact that De = eD = Iwe also have 



(x) = ^"p'ji^) = ^[Yl ^^^^^ ) 

\fc>0 / 



k>0 



Combining the above results, we have 

N-l N-1 
j=0 j=0 k>0 

= "'k,2j+i¥^k (x) e ip2j+i (x) 

j=0 k>0 
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Note that afc,2j+i = unless \k — {2j + 1)| = 1, that is unless k is even. Thus we can 
rewrite the sum as 



N-l 



^2JV,2JV+l 



j=0 



k,j>0 
k even 
k<2N 



i>o 



A:>0 
A: even 
k<2N 



where the last term takes care of the fact that we are counting an extra term in the 
sum that was not present before. The sum over j on the right is just (pk{y), and 



a2^,2iv+i = -\/^- Therefore 



Af-l 



Y '^2i+i i^) e v^2j+i (y) = 5^ v'fc (x) ^k (y) 

j=0 



k>0 
k even 
k<2N 

N 



2N + 1 



^2N {X)^^2N+1 (y) 



Y ^'ij (^) ^'^j (y) - \l —w— V2N {x) e ¥^2^+1 {y) 

3=0 



It follows that 



■ 2N 

.j=o 



2N + 1 



We redefine 



2N 



Sn{x, y) = Y'^r^ (x) iv) = SnIv, x) 



(4.2.8) 



n=0 



SO that the top left entry of Kj^{x,y) is 



'2N + 1 

Sn{x, y) + \/ — - — (x) e ip^^^^ {y) . 
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SnDJ = I S{x,y)-^f{y)dy 



■^SN{x,y) ) f{y)dy, 



so that —j;^SN{x,y) is in fact the kernel of SnD- Therefore ()4.1.4j) now holds with 
being the integral operator with matrix kernel K^^^j^^x^y) whose (i,j)-entry 
^Jjv^a;,?/) is given by 



Ki^Nix,y) 

K^Ni^^y) 



SN{x,y) + 
SDN{x,y) 
SNix,y) + 
SN{x,y) + 



2iV + 


1 







d ( 




dy y 




2iV + 


1 


2 


2N ^ 


1 


2 



^2N {X)^^2N+1 (y) 



2N+1 



2 r 2N K'^ J r 2]V+1 

f2N ix)^f2N+l (y) 



We let 2N+ 1 ^ N so that is assumed to be odd from now on (this will not matter 
in the end since we will take 00). Therefore the K^^'^{x, y) are given by 



Kf^Ni^^y) 

K4,N\x,y) 



SN{x,y) + \l — ip^_^{x)eip^{y) 



SDNix,y) - \ —ip^_^{x)(p^{y) 



SN{x,y) + \l — ip^_^{x)eip^{y) 



SN{x,y) + \l — eip^{x) ip^_,iy) 



where 



N-l 



SN{x,y) = ^ip^{x)ip^{y). 



n=0 
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Define 



so that 



Notice that 



Therefore 



1/4 



N 



1/4 



i,Ni^^y) = 7;^i^) [SN{x,y) +ip{x)eif{y)] X{y), 



K. 



Kt:M\x,y) = \x{x) [SDN{x,y)-ij{x)^{y)] X{y), 

.(2,1). . 1 



K 



(2,2), 
4,JV ' 



{x.y) = - X{x) [eSNix,y) + eip{x)eip{y)] X{y), 
x^y) = \ [SN{x,y) + e^{x)%l){y)] X{y). 



■X{S + ^® e^) X = K^^;i\x,y), 



X {SD -i)® ^)X 



K'i'i\x,y), 



^X {eS + e^® eip) X = K^^i\x,y), 



-X {S + eip® etp) X 



KfPix.y). 



4,JV 



X 



^ S + i)® eip SD - tp ® ip^ 



V 



X. 



I 



(4.2.9) 



eS + eip®eip S + eip®ip 

Note that this is identical to the corresponding operator for /? = 4 obtained by Tracy 
and Widom in |23|; the only difference being that ip, and hence also S*, are redefined 
to depend on A. This will affect boundary conditions for the differential equations we 
will obtain later. 



4.3 Edge scaling 

4.3.1 Reduction of the determinant 

We want to compute the Fredholm determinant (j4.1.4p with K^^j^ given by (|4.2.9|) 
and / = X(t,oo)- This is the determinant of an operator on L'^{J) x L'^{J). Our first 
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task will be to rewrite the determinant as that of an operator on L^{J). This part 
follows exactly the proof in j^Sj • To begin, note that 

[S,D]=ip®^P + ^P®(p (4.3.1) 

so that, using the fact that D e = e D = I , 

, S] = eS — S e 

= eSDe-eDSe = e [S,D] e 
= eip ® ifj e + eip ® if e 
= e V5 ® eV + e ® 

= -eip ® e^p - e^p IS) e(f, (4.3.2) 
where the last equality follows from the fact that e* = — e. We thus have 

D {eS D - eij®ip) = S D - ^ ip. 

The expressions on the right side are the top matrix entries in ()4.2.9|) . Thus the first 
row of -ft'4,Ar is, as a vector, 

D {eS + t^p ® eip.eS D - e^p ® Lf) . 

Now PA2|1 implies that 



e S + eip ® e^p = S e — e^p ® eip. 
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Similarly ()4.3.H) gives 



e [S ,D] = e^p i^ip + eijj ® if, 



so that 



eSD-e'ip(g)(p = eDS + e(p(g)'ilj = S + e(p(S)ip. 



Using these expressions we can rewrite the first row of K4,^n as 



D {Se-eip®eip,S + eip®^) 



Now use ()4.3.2|) to show the second row of K/^^jq is 



Therefore, 



X 
/ 



\ 



XD 
X 



S e — eip ® eip 
\ 



X 



J 



J 



S + e(p<^tp 

{Se- eLp®e'(p) X {S + eLp®i))X 
{Se- eip®ei))X {S + eip®i})X 



Since A'4,Ar is of the form AB, we can use 12. 1.2] and deduce that _D4 Ar(s, A) is un- 
changed if instead we take ^ to be 



{Se-eip®ei/j) X {S + e(p®i/j)X 
y{Se-eip®eij) X {S + eip®ij) X J y 

^ {Se-eip®ei/j) XD {S + e(p®-^)X^ 



^ ^ XD ^ 



X 



V 



{S e - e^p etp)XD {S + eip®i))X 



J 
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Therefore 

\ 



-D4,Ar(s, A) = det 



I - ^ {Se- e(p®eij) \XD {S + e^p ® ip) XX 

-\ {Se- e^^e^j) \XD I - ^ {S + etp (g) ip) XX ^ 



(4.3.3) 

Now we perform row and column operations on the matrix to simphfy it, which do 
not change the Fredholm determinant. Justification of these operations is given in 
[22] • We start by subtracting row 1 from row 2 to get 

( I - \ {Se-e^®eil)) \XD -\ {S + ® i)) \x\ 
-I I 



Next, adding column 2 to column 1 yields 

( I -\{Se-e^®ei)) \XD-\{S + e^®i}) \X -\ {S + ® i)) \x\ 

/ 



Thus the determinant we want equals the determinant of 



I- - {Se-e^®e^lj) XXD- - {S + eLp®ij) XX (4.3.4) 
2 2 



So we have reduced the problem from the computation of the Fredholm determinant 
of an operator on L'^[J) x L'^[J), to that of an operator on L'^{J). 

4.3.2 Differential equations 

Next we want to write the operator in ()4.3.4j) in the form 

(/ - K2,n) Pij , (4-3.5) 
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where the and f3i are functions in L'^{J). In other words, we want to rewrite the 
determinant for the GSE finite dimensional perturbation of the correspond- 

ing GUE determinant. The Fredholm determinant of the product is then the product 
of the determinants. The hmiting form for the GUE part is aheady known, and we 
can just focus on finding a hmiting form for the determinant of the finite dimen- 
sional piece. It is here that the proof must be modified from that in j2S]- A little 
rearrangement of ()4.3.4j) yields (recall e* = — e) 



I-^SX-^SeXD-^ecp (g) Xip-^ecp ® ipeXD. 



Writing e [X , D] + X for eX D and simplifying gives 



I - \SX- XeLf ® i)X- - Se [X,D] --eip ® i) e [X,D'_ 

2 2 



Let y/X^p — >■ V?, and y/Xip ^ ijj so that \ S ^ S and ()4.3.4|) goes to 



I - SX-eip ® tpX-^Se [X,D]-^eif (g) ip e [X,D'_ 



Now we define R:= [I — S X)~^ SX = (I — S XY^ — I (the resolvent operator of SX), 
whose kernel we denote by R{x,y), and := (/ — SX)"^eip. Then fj4.3.4p factors 
into 

A={I-SX)B. 

where B is 



I-Q,0Xi;-^{I + R)Se[X,D]-^{Q,®i;)e[X,D' 



Hence 



D4,7v(t,A) =D2,;v(t,A) det(5). 
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In order to find det(-B) we use the identity 



2m 



6 =5^(-l)'=efc®5fc, (4.3.6) 



k=l 



where and Sk are the functions e{x — a^) and 5{x — a^) respectively, and the are 
the endpoints of the (disjoint) intervals considered, J = U^^^(a2k-i, 0.2 k)- In our case 
m = 1 and ai = t, a2 = 00. We also make use of the fact that 

aiS)b-c®d={b,c)-a®d (4.3.7) 

where ( . , . ) is the usual L^-inner product. Therefore 

2 

{Q,®^)e[X,D] = J2i-'^)'Qe^i^-ek^Sk 

k=l 
2 

= J](-l)'=(^,efc) Q,® 



fe=i 

It follows that 

A.jv(t,A) 
/^2,iv(t,A) 

is the determinant of 



(4.3.8) 



1 2 

I-Q,®X^--J2i-'^)'liS + RS)ek + {tlJ,ek) Qe] ® 4- (4.3.9) 



2 

k=l 



We now specialize to the case of one interval J = {t,oo), so m = 1, ai = t and 
a2 = 00. We write ej = ei, and €00 = €2, and similarly for 6k- Writing out the terms 
in the summation and using the fact that 



eoo = (4.3.10) 
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yields 



I-Q,®Xij + ^ [iS + RS)et + {^lj,et) Qe]®St + ^ [{S + RS) 1 + {-iP ,1) Q,] ® 5oo 

(4.3.11) 

Now we can use the formula 

det (^I-Y,ai® A j = det {6jk -(«,-, pk ))i<j,k<L (4.3.12) 

In order to simplify the notation in preparation for the computation of the various 
inner products, define 



Q{x,X,t) := {I-SX)-^ip, P(x,A,t) := (I-SX)-^^, 

Q,{x,X,t) := {I - Sxy^if, P,{x,X,t) := {I - SX)-\iIj, (4.3.13) 



tie : = 



Q{t,\,t), 

Qeit,X,t), 

{Q,Xeip) = {Q,,Xip) 
{P,Xeip) = {Qe,Xip). 



Pn ■-- 

Pe ■-- 



P{t,X,t), 

Pe{t,X,t), 

(P,Xe^) 



Pe,X^) 



(4.3.14) 



Pi '■= / et{x) P{x,t) dx, Q4 := / et{x) Q{x,t) dx, 71^:= / et{x) R{x,t) dx, 
Jr Jr Jr 

(4.3.15) 

where we remind the reader that et stands for the function e(x — t). Note that all 
quantities in ()4.3.14p and ()4.3.15p are functions of t and A alone. Furthermore, let 

= e(p{oo) = - / (p{x)dx, = eij{oo) = - / ij{x)dx. (4.3.16) 
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Recall from the previous section that when /3 = 4 we take to be odd. It follows 
that (f and ip are odd and even functions respectively. Thus when f] = 4, = 
while computation using known integrals for the Hermite polynomials gives 

c, = (vriV)V^2-3/-^/^i^v^. (4.3.17) 



limc^ = W^. (4.3.18) 



Hence computation yields 

Af^'i) "'^ V 2 

At t = oo, 

Ue(oo) = 0, Q'e(oo) = (4.3.19) 

^4(00) = -c^, Q4(oo) = -c^, 7^4(oo) = 0. (4.3.20) 
In (14.3. 12|1 . L = 3 and if we denote 04 = ( t/^ , ), then we have explicitly 

«i = Q., «2 = -^ [{S + RS)et + a^Qe], «3 = [{S + RS)1 + ,1) Q,] , 
However notice that 

{{S + RS) et,6^) = {et,6o^) = 0, {{S + R S) 1 ,6^) = {1 , R^) = (4.3.21) 

and ( , ^oo ) = = 0. Therefore the terms involving = S^o are all and we can 
discard them reducing our computation to that of a 2 x 2 determinant instead with 



«i = Qe, «2 = -- [{S + RS)et + a, Q,] , A = X^j, P2 = (4.3.22) 
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Hence 



>1 ) = fe, (ai ,/32 ) = Qe, 

! "2 , A ) = (^4 - ^4 +04 Se) 

>2,/32) = (7^4 + a4g,). 



(4.3.23) 
(4.3.24) 
(4.3.25) 



We want the limit of the determinant 



det {5jk - {aj,pk))i<j. 



k<2 ' 



(4.3.26) 



as ^ oo. In order to get our hands on the hmits of the individual terms involved 
in the determinant, we will find differential equations for them first as in [22] • Adding 
04/2 times row 1 to row 2 shows that 04 falls out of the determinant, so we will not 
need to find differential equations for it. Thus our determinant is now 



det 



/ 



\ 



\ 



\Va l + |7^4 



(4.3.27) 



Proceeding as in j2S] we find the following differential equations 



d_ 

Tt 

d ^ 



-q^ (7^4 + 1) , 



d 

Tt 

j-p, = -p, in, + 1] 



(4.3.28) 
(4.3.29) 
(4.3.30) 



Now we change variable from t to s where t = r(s) = V2N + s/{y/2 N-^^^) and take 
the limit N ^ 00, denoting the limits of q^, 7^4, Q4, TZ,, and the common limit of 
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Me and respectively by q, V^, Q4, 7^4 and u. Also P4 and Q4 differ by a constant, 
namely ^4 = ^4 + V2/2. These limits hold uniformly for bounded s so we can 
interchange limAr^oo and Also liraN^oo N'^^^q^^ = limAr^oo A^~^^V]v = 1 5 where 
q is as in (jl.l.lOj) . We obtain the systems 



iLp, = __L,(K, + l). iL,^,__L,(2p,H.y|j, (4.3.32) 

The change of variables s ^ fi = q{x, X) dx transforms these systems into constant 
coefficient ordinary differential equations 



-^u =^q, -^g = (1 - 2 n) , (4.3.33) 



Since lim^^oo A* = 0, corresponding to the boundary values at t = 00 which we found 
earlier for V4,TZa, we now have initial values at yU = 0. Therefore 

^(/i = 0) = g(/i = 0) = 0, (4.3.35) 



V^il^ = 0) 



7^4(/x = 0) = 0. 



(4.3.36) 



CHAPTER 4. THE M^" LARGEST EIGENVALUE IN THE GSE 
We use this to solve the systems and get 



49 



7^4 



u = - 



q = 




(4.3.39) 



(4.3.38) 



(4.3.40) 



(4.3.37) 



Substituting these expressions into the determinant gives ()1.2.2j) . namely 



D4(s, A) = D2{s,X) cosh^ 



/^(g, A) 
2 



) 



(4.3.41) 



where Djs = limAr^oo -D/j^at. Note that even though there are A-terms in ()4.3.39p and 
fl4.3.4()|l . these do not appear in the final result ()4.3.4H1 . making it similar to the GUE 
case where the main conceptual difference between the m = 1 (largest eigenvalue) 
case and the general m is the dependence of the function g on A. The right hand 
side of the above formula clearly reduces to the (3 = 4 Tracy- Widom distribution 
when we set A = 1. Note that where we have D^^s^X) above, Tracy and Widom 
(and hence many RMT references) write ^4(5/1/2, A) instead. Tracy and Widom 
applied the change of variable s sj \/2 in their derivation in [221 agree with 

Mehta's form of the /? = 4 joint eigenvalue density (see ()A.1.3j) and footnote in the 
Introduction), which has — 2x^ in the exponential in the weight function, instead of 
—x^ in our case. To switch back to the other convention, one just needs to substitute 
in the argument s/\/2 for s everywhere in our results. At this point this is just a 
cosmetic discrepancy, and it does not change anything in our derivations since all the 
differentiations are done with respect to A anyway. It does change conventions for 
rescaling data while doing numerical work though (see Section rA.4jl . 



Chapter 5 

The distribution of the m^^ largest 
eigenvalue in the GOE 



5.1 The distribution function as a Fredholm deter- 
minant 

The GOE corresponds case corresponds to the speciahzation /? = 1 in so that 



Gi,^(t,A) = cf^ / ■■■ / ni^:^ -^^l + fi^j)) dxi---dxN (5.1.1) 



where w{x) = exp (— x^), f{x) = — AXj(x), and Cj depends only on A^. As in the 
GSE case, we will lump into C['^^ any constants depending only on that appear in 
the derivation. A simple argument at the end will show that the final constant is 1. 
These calculations more or less faithfully follow and expand on [25- We want to use 
()2.1.2|) . which requires an ordered space. Note that the above integrand is symmetric 
under permutations, so the integral is n\ times the same integral over ordered pairs 




TV 



N 
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Xi < . . . < xat. So we can rewrite 15.1.11 as 



„ „ N N 

(iV!) /■■■/ Yl{xk-Xj)Y[w{xk)Yl{l + f{xk))dxi---dxN, 



where we can remove the absolute values since the ordering insures that {xj — Xj) > 
for i < j. Recall that the Vandermonde determinant is 

An{x) = det(xr^)i<^,,<^ = (-1)^^^ Y[{xj - Xk). 

j<k 

Therefore what we have inside the integrand above is, up to sign 

det(xr^M;(xfc) (1 + /(xfc))),<^, 



Note that the sign depends only on A^. Now we can use (j2.1.2|) with (pj\ 



X] 



x^ ^ wix) 



[1 + f{x)). In using ()2.1.2p we square both sides so that the right hand 



l<j,k<N 



side is now a determinant instead of a Pfaffian. Therefore G\^{t, A) equals 

Shifting indices, we can write it as 

C[^^det( sgn{x -y)x^ y''{l + f{x))w{x)w{y)dxdy) (5.1.2) 



l<j,fc<JV-l 



where cl^"* is a constant depending only on A^, and is such that the right side is 1 
if / = 0. Indeed this would correspond to the probability that Xp'^^^'^^ < oo, or 
equivalently to the case where the excluded set J is empty. We can replace x^ and y'' 
by any arbitrary polynomials Pj{x) and Pk{x), of degree j and k respectively, which 
are obtained by row operations on the matrix. Indeed such operations would not 
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change the determinant. We also replace sgn(x — y) by e{x — y) = | sgn{x — y) which 
just produces a factor of 2 that we absorb in C^'^K Thus G\j^{t, A) now equals 



0<j,fc<JV-l 

" (5.1.3) 



Let ipj{x) = Pj{x) w{x) so the above integral becomes 

det( f [ e{x-y)ij,{x)My){'^ + f{^) + f{y) + f{^)f{y))dxdy 



0<j,k<N-l 

" (5.1.4) 



Partially multiplying out the term we obtain 




+ / i ^{x-y)i^j{x)^k{y){f{x) + f{y) + f{x)f{y))dxdy 



0<j.k<N-l 

Define 

f e{x-y)^tijj{x)My)dxdy] , (5.1.5) 

so that Gf ^(t, A) is now 

det (m+ [ [ e{x-y)^,{x)My)ifi^) + fiy) + f{x)f{y))dxdy 

' 0<j.k<N-l 

Let e be the operator defined in ()4.1.2|) . We can use operator notation to simplify 
the expression for Gl^(t,X) a great deal by rewriting the double integrals as single 
integrals. Indeed 




e{x - y)i/jj{x)i/jk{y) f{x)dxdy = / f{x)ipj{x) J e{x - y)i)k{y) dy dx 

fi)jei)kdx. 
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Similarly, 





e{x - y)ijj{x)iJk{y) fiy)dxdy = - / j e{y - x) iljj{x) ipki^y) f{y) dx dy 

f{y)i^kiy) / e{y - x)ipj{x)dxdy 



f{x)iJk{x) J e{x-y)^ljj{y)dydx 
ftpketpidx. 



Finally, 





{x - y) i)j{x) i)k{y) f{x) f{y)dxdy 

e{y - x) ipj{x) ipk{y) f{x) f{y)dxdy 
f{y)i^k{y) j e{y - x) f{x)^jjj{x) dxdy (5.1.6) 
f{x)tpkix) J e{x -y) f{y)i)j{y)dydx 



It follows that 



Gl^it,X)=Cf^ det(M+ / [f^^eil;k-fipke^j-fipke{f^j)]d 



X 



0<j,fc<JV-l 

- -(5.1.7) 



If we let M ^ = if^jk)^^ . J, and factor det(M) out, then Gl j^(t, A) equals 
Cf^ det(M) det 

M-'-( [ [f^,e4^k-f4^keiJj-fiJke{f^,)]dx) ] (5.1.8 



0<j,k<N-l , 



0<j,k<N-l 



where the dot denotes matrix multiplication of M ^ and the matrix with the integral 
as its (j, /c)-entry. define rjj = J2j f^j k ipk and use it to simplify the result of carrying 
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out the matrix multiplication. From ()5.1.5|) it follows that det(M) depends only on 
TV we lump it into C[^\ Thus Glj^{t, A) equals 



Cj^^detl/+( / [fr]jeiJk-fipker]j-f'ilJke{frij)]dx 



0<j,k<N~l , 



0<j.k<N-l 

(5.1.9) 



Recall our remark at the very beginning of the section that if / = then the integral 
we started with evaluates to 1 so that 



Cl""^ det(/) =Cl 



(TV) 



(5.1.10) 



which implies that c[^'^ = 1. Now G^^(t, A) is of the form det(J + AB) where 
A : L2(J) X L2(J) ^ is a iV X 2 matrix 



whose j^^ row is given by 



A 



A, 



= Ajix) = (-/ er]j- fe if rjj) f rjj) . 



( \ 



Therefore, if 

/ 

9 = 

then Ag is a, column vector whose j^^ row is {Aj ,g 



e L\J) X L\J), 



iA9)j= / [-f^Vj-f^ifVj)] 9idx+ / fr]jg2dx 
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Similarly, B : ^ L'^{J) x L^{J) is a 2 x N matrix 



B={ B, B2 ... B 



N ) ' 



whose j*^ column is given by 



Bj = Bj{x) 



Thus if 



e C 



N 



then B h is the column vector of L'^{J) x L'^{J) given by 



Bh 



Clearly A B : ^ and B A : L'^{J) x L^J) L'^{J) x L\J) with kernel 



Hence / + BA has kernel 
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which can be written as 

Since we are taking the determinant of this operator expression, and the determinant 
of the second term is just 1, we can drop it. Therefore 





( I 


°1 


/ 







( 



Glj^{t,\) = det 



det(/+ i^i,iv/), 



N 



where 

y - Ei e^Ai er]j-e Ej^^j ® Vj j 
and -ft'i^AT has matrix kernel 



We define 



- Ei,fc t^jkeMy) T.j,k ^A^) f^jk My) 
y - J2j,k e i^jix) Hk e My) - - y) J2j,k e i^ji^) ^^jk My) J 



N-1 



SN{x,y) = - ^ i)j{x)^ijkeMy)- 

j,k=0 

Since M is antisymmetric, 

Af-l N-1 N-1 

SN{y,x) = - ^ ipj{y)iijk^i)k{.x) = ijj{y)fikjeiJkix) = eipj{x)iijk'ilJk{.y)- 

j,k=0 j,k=Q j,k=0 
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Note that 



N-l 



eSNix,y) = ^ eipj{x) fijkeipkiy), 

j,k=Q 



whereas 



N-l 



y j,k=o 



So we can now write succinctly 



Ki,N{x,y) 



^ SN{x,y) -^SN{x,y) ^ 



V 



eSN{x,y)-e SN{y,x) 



J 



So we have shown that 



where 



(5.1.11) 



(5.1.12) 



Gi,N{t,X) = y^/^i,^ (t,A) 



Di,^(t,A) =det {I + K^^Nf) 



where A'i^tv is the integral operator with matrix kernel Ki^fyf{x, y) given in ()5.1.11|) . 



5.2 Gaussian specialization 



We specialize the results above to the case of a Gaussian weight function 



w{x) = exp (— 



(5.2.1) 



and indicator function 



f{x) = -\X, J=(t,oo) 



Note that this does not agree with the weight function in (jl.l.lll . However it is a 
necessary choice if we want the technical convenience of working with exactly the same 
orthogonal polynomials (the Hermite functions) as in the /3 = 2,4 cases. In turn the 
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Painleve function in the limiting distribution will be unchanged. The discrepancy is 
resolved by the choice of standard deviation. Namely here the standard deviation on 
the diagonal matrix elements is taken to be 1, corresponding to the weight function 
()5.2.H) . In the /3 = 2, 4 cases the standard deviation on the diagonal matrix elements 
is 1 / \/2, giving the weight function ()4.2.1|) . We expand on this in Section I A. 41 Now 
we again want the matrix 



M 




e{x - y) %l)j{x) 'ijjk{y)dxdy 



0<j,k<N-l 



%l)j{x)ei)k{x)dx 



0<j,k<N-l 



to be the direct sum of ^ copies of 



^ l^ 



v 



-1 



so that the formulas are the simplest possible, since then fijk can only be or ±1. In 
that case M would be skew-symmetric so that = —M. In terms of the integrals 
defining the entries of M this means that we would like to have 



and otherwise 



'4'2m+l{x) eip2nix)dx = -5m,n, 



ipj(x) ipkix) dx = 0. 
dx 



It is easier to treat this last case if we replace it with three non-exclusive conditions 



'Ip2m{x)eip2nix)dx = 0, 
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(so when the parity is the same for j, k, which takes care of diagonal entries, amon^ 
others), and 

^j{x) eipkix) dx = 0. 



whenever \j — k\ > 1, which targets entries outside of the tridiagonal. Define 



1 



(Pn{x) = —Hn{x) exp{-x /2) for Cn = \J '2Pn\^/n 



where the if„ are the usual Hermite polynomials defined by the orthogonality condi- 
tion 



It follows that 



Now let 



j Hj{x) Hk{x) e dx = c] 6j,k- 



(Pj{x) (fk{x) dx = 6j^k- 



V'.^^, (x) = ^ ( x) (x) = ( x) . (5.2.2) 



This definition satisfies our earlier requirement that ipj = pj w for 

w{x) = exp (— x^/2) . 

In this case for example 

P2n{x) = —H2n (x) . 

CtJ, 
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With e defined as in ()4.1.2j) . and recalling that, if D denote the operator that acts by 
differentiation with respect to x, then De = eD = I, it follows that 

/d 
'P2m. {x)e — ip,^^, {x)dx 

'^2m (X) <^2n + l i^) d {x) 



5. 



m,ni 



as desired. Similarly, integration by parts gives 



i^2ra+Axyi^2n{x)dX= / ^ {x) € ^ {x) d X 

¥^2^ (a;) V2r. {x)dx 

- / <^2,„ ^2n + l d (X) 



Also is even since and (f^n are. Similarly, ip2n+i is odd. It follows that eip^n, 
and e^Aan+i' respectively odd and even functions. From these observations, we 
obtain 

^2A^)^^2Mdx = 0, 



since the integrand is a product of an odd and an even function. Similarly 

A^+li^)^'^2n. + li^)dX = 0. 

Finally it is easy to see that if \j — k\ > 1, then 

■ipj{x) eipf,{x)dx = 0. 
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Indeed both differentiation and the action of e can only "shift" the indices by 1. Thus 
by orthogonahty of the (fj, this integral will always be 0. Thus by our choice in 
()5.2.2p . we force the matrix 



M 



0<j,k<N-l 



to be the direct sum of ^ copies of 



^0 1^ 



This means M ^ = —M where M ^ = {f-tj,k}- Moreover, /ij^^ = if j,k have the 

same parity or \j - A;| > 1, and fJ,2j,2k+i = ^jk = -fJ'2k+i,2j for j, /c = , ^ - 1- 

Therefore 



N-l 

SN{x,y) = - ^ ipj{x)njkeilJkiy) 

j,k=Q 

N/2~l N/2-1 

= - ip2j{x)eilj2j+i{y) + i/j2j+i{x)eip2j{y) 

j=0 j=0 
'n/2-1 N/2~l ^ 

j=0 j=0 



Manipulations similar to those in the /9 = 4 case (see ()4.2.3|1 through ()4.2.8j) ) yield 



SN{x,y) 



'N-l 



.j=o 



We redefine 



N-l 



Sn{x, y) = Y = 
i=o 
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so that the top left entry of Ki^]\f{x,y) is 



Swix, y) + \l^ Vn^i i.^) (e Vn) (y) 



If Sn is the operator with kernel Sn{x, y) then integration by parts gives 



SNDf = J Six,y)-^fiy)dy 



■^SN{x,y) ] f{y)dy, 



so that —j^SN{x,y) is in fact the kernel of SnD. Therefore (j5.1.12p now holds with 
i^i^TV being the integral operator with matrix kernel Ki^i<^{x,y) whose (2,j)-entry 
K[''^\x,y) is given by 



Define 



so that 



K[]i\x,y) = 

K[%\x,y) = e 
K[f{x,y) = 



Swix, y) + \l — (^) (e Vn) (y) 



SDN{x,y) - ^ (x) (eip^) (y) 



Sn{x, y) + \l — (^) (e Vn) (y) - 1 



SN{x,y) + V y (eV'^) (x) (y) 



1/4 

V^{x), ip{x) = ( y 1 ^iv-i(^)' 



K[];i\x,y) = Xix) [SN{x,y)+^{x)eipiy)] X{y), 
K[]^\x,y) = Xix) [SD^ix,y)-^ix)^iy)] 

{^,y) = X{x) [eSN{x,y) + e^j^x) eip{y) - e{x - y)] X{y), 
Kf^\x,y) = X{x) [SN{x,y) + eip{x)^{y)] X{y). 
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Note that 



Hence 



X (S + ih ® 


eip) X = 


1,N 


'x. v). 


X{SD-^^ 


^ ^)x = 


-"■1,7V 




X {eS + eip ® e^p 


-e)X = 


-"■1,7V 






e^) X = 


-'^^1,7V 





( 



N 



X 



\ 



\ 



X. 



J 
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S + tpiS) eip SD - if 
e S + eip ® e(p — e S + eip®eip 

Note that this is identical to the corresponding operator for (3 = 1 obtained by 
Tracy and Widom in [22], the only difference being that (p, ip, and hence also 5", are 
redefined to depend on A. 



5.3 Edge scaling 

5.3.1 Reduction of the determinant 

The above determinant is that of an operator on L'^{J) x L'^{J). Our first task will 
be to rewrite these determinants as those of operators on L'^{J). This part follows 
exactly the proof in jlH] . To begin, note that 



[S ,D] =ip®ilj + ilj®ip 



(5.3.1) 
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[6,5*] = eS — S e 

= eSDe-eDSe = e [S,D] e 
= eip ® ifj e + eip ® if e 

= —e(p®e^lj — etjj'S>e(p^ (5.3.2) 
where the last equahty follows from the fact that e* = — e. We thus have 



D {eS + etp (g) eip) = S + i)®eip, 
D {eSD -etp0ip) = S D - ^ ip. 

The expressions on the right side are the top entries of Ki n. Thus the first row of 
Ki n is, as a vector, 

D {eS + ei/j e(p,eS D - eijj if) . 
Now ()5.3.2j) implies that 

eS + eip^eip = Se — eip^eip. 

Similarly ()5.3.1|) gives 

e [S ,D] = e(p ®ip + eijj ® (f, 

so that 

eSD-eip0Lp = eDS + eLp0'ilj = S + eLpiS)^. 
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Using these expressions we can rewrite the first row of Ki^n as 



Applying e to this expression shows the second row of Ki^n is given by 



{eS -e + ei)®eip,S + eip®i)) 



Now use ()5.H.2|) to show the second row of Ki^]^ is 



Therefore, 



N 



X 

/ 



\ 



Se — e + eip^eip S + e^p®ilj 
\ ( 



X 



J 



XD 
X 



J 



\ 



{Se-eif^eip) X {S + e(p®ilj)X 
y {Se-e + e<^(^etp)X {S + e^(^ij)X J 



Since Ki^^ is of the form AB, we can use the fact that det(/ — AB) = det(/ 
and deduce that Di i\f{s, A) is unchanged if instead we take Ki^^ to be 



K 



l.N 



f 



{Se- eif®ei)) X {S + eip®i))X 
{Se-e + eLp®eilj)X {S + eLp®ilj)X 



^ ^ XD ^ 



V 



X 



^ {Se-e^®ei))XD {S + e^®i))X^ 



{Se-e + eip®e^)XD {S + eip®il))X 



J 



CHAPTER 5. THE M^" LARGEST EIGENVALUE IN THE GOE 



66 



Therefore 



Di^n{s, A) = det 



^ I - {Se- eip(S)eip) XXD - {S + eif (g) ip) XX ^ 



v 



{Se-e + e(p®etp) XXD I - {S + e^p <S> ij) XX 



. (5.3.3) 



Now we perform row and column operations on the matrix to simphfy it, which do 
not change the Fredholm determinant. Justification of these operations is given in 
[22] • We start by subtracting row 1 from row 2 to get 



I - {S e - eip etfj) XXD - [S + e^p ® tp) XX 
-I + eXXD I 



Next, adding column 2 to column 1 yields 



^ I -{Se-eif^e-^) XX D - {S + eip ^j) XX -{S + eip0ij) XX ^ 



V 



eXXD 



I 



Then right-multiply column 2 by —eXXD and add it to column 1, and multiply row 

2hyS + eip®il) and add it to row 1 to arrive at 



det 



^ I - {Se-e^®eil)) XXD + {S + e^®il)) XX {eXXD - I) 0^ 



v 







Thus the determinant we want equals the determinant of 



I -{Se-e^®ei)) XXD + {S + e^®^) XX{eXXD-I). 



(5.3.4) 



So we have reduced the problem from the computation of the Fredholm determinant 
of an operator on L'^{J) x L'^{J), to that of an operator on L'^{J). 
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5.3.2 Differential equations 

Next we want to write the operator in ()5.3.4j) in the form 

(/ - K2,n) - 5^ ® A j . (5-3.5) 

where the Oj and /3j are functions in L'^{J). In other words, we want to rewrite the 
determinant for the GOE finite dimensional perturbation of the correspond- 

ing GUE determinant. The Fredholm determinant of the product is then the product 
of the determinants. The hmiting form for the GUE part is aheady known, and we 
can just focus on finding a hmiting form for the determinant of the finite dimen- 
sional piece. It is here that the proof must be modified from that in A little 
simplification of ()5.H.4|1 yields 

I - \SX-XS {I - \X) eXD -\ {e^ ® Xi))- \ {e^ ® - \X) eXD. 

Writing e [X , D] + X ioi eX D and simplifying (1 — A X) X to (1 — A) X gives 



I - \SX - \ {I - X) SX - \ {eif ® X^j) - \ {I - \) {eif ® X^j) 
- \ S {I - \ X) e [X , D] - \ {e^ ® ^l)) [l ~ \ X) e [X , D] 

= I - {2\ ~ \^) SX - (2\ ~ \^) {e^ ® X^^) - \ S {I - \ X) e [X , D] 
-\{eip ® ^){l-\X)e [X,D]. 
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Define A = 2 A — and let if ^ if, and ^/Xijj ^ ijj SO that A 5 ^ 5 and (pnOjl 
goes to 

I-SX-{eip ® Xij) -ts{l - XX)e[X,D] 

A 

- i (ev9 (g) ip) (1 -XX)e [X,D]. 
A 

Now we define R := {I — S X)~^ S X = {I — S X)~^ — I (the resolvent operator of S X), 
whose kernel we denote by R{x,y), and := {I — SX)~^eip. Then fj5.3.4p factors 
into 

A={I-SX)B. 

where B is 



I-{Q, ® Xip) - ^ {I + R) S {I - \ X) e [X , D] 
A 

-jiQe^ ^){l-XX)e[X,D], A^l. 
A 

Hence 

Di,^(s,A) =D2,;v(s,A) det(5). 

Note that because of the change of variable A S* ^ S*, we are in effect factoring 
/ — (2A — A^) S, rather that J — A S* as we did in the /3 = 4 case. The fact that we 
factored / — (2A — A^) S* X as opposed to / — A S* X is crucial here for it is what makes 
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B finite rank. If we had factored I — XS X instead, B would have been 

2 

B = I-X + {I - XX)ek^6k- X {I + R) ® X^/^ 

k=l 
2 

-A 5^(-l)'=(^,(/-AX) efc) {{I + R) e^)®6k 

k=l 

- X{1 - X) {S + RS) X - X{1 - X) {{I + R) e^) ®X^ 

The first term on the last line is not finite rank, and the methods we have used 
previously in the /3 = 4 case would not work here. It is also interesting to note 
that these complications disappear when we are dealing with the case of the largest 
eigenvalue; then is no differentiation with respect to A, and we just set A = 1 in all 
these formulae. All the new troublesome terms vanish! 
In order to find det(i?) we use the identity 

2m 

e =5^(-l)'efc®4, (5.3.6) 

k=l 

where and 6k are the functions e(x — ak) and 6{x — au) respectively, and the ak are 
the endpoints of the (disjoint) intervals considered, J = U^=i{a2k-i,CL2k)- We also 
make use of the fact that 



a(Sib-c^d= {b,c) -a^d (5.3.7) 



where ( . , . ) is the usual L^-inner product. Therefore 



2m 

(g,®^)(l-AX)e = ^(-l)'=g,®V^- (1- AX)efc®4 

k=l 
2m 

= J](-l)'(^,(l-AX)efc) Q,® 5k. 
fe=i 
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Di,n{s, A) 
D2,n{s, A) 



equals the determinant of 



, 2m 

I-Q,m^-~ [{S + RS) {1 - XX)ek + ,{1 - XX)ek) Q,] ®4. (5.3.8) 

^ k=i 

We now specialize to the case of one interval J = (t, oo), so m = 1, ai = t and 
a2 = oo. We write et = €i, and 600 = ^2? and similarly for 6k- Writing the terms in 
the summation and using the facts that 



Coo = (5.3.9) 



and 



then yields 



(1 - A X) et = -i (1 - A X) + (1 - A X) X, (5.3.10) 



I-Q,®Xij-^ [{S + RS) {1 - XX) + (iP , {1 - XX)) Q,] (g) i6t - 6^) 
2A 

+^[iS + RS)il-XX)X + iij,il-XX)X) Q,](g)6t 
X 



which, to simplify notation, we write as 



I - Q, ^ X^P - ^ [{S + RS) {1 - XX) + ai,A Qe] ® (St - 5oo) 
2X 



+^ [{S + RS){1-XX)X + ai,A Qe] ® St, 
X 



where 

ai,A = ( ^ , (1 - A X) ) , fii,A = ( ^ , (1 - A X) X ) . (5.3.11) 
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Now we can use the formula: 

det ^/ - ^ tti (g) = det {5jk - {aj , (3k ))i<j^k<L 
In this case, L = 3, and 

ai = Qe, «2 = Y [{S + RS){1- XX) + ai,xQe], 
A 

as = -^[{S + RS)il-XX)X + ai,A Qe] , 
A 

/5i = Xip, ^2 = St- Soo, (3-3 = St.. 
In order to simplify the notation, define 

Q{x,X,t) := (I-SX)-^^, P{x,X,t) := {I-SX)-^ij, 

Qe{x,X,t) := {I -Sxy\ip, Pe{x,X,t) := {I-SX)~\^Ij, 



:= Qeit,X,t), 

Ue := {Q ,Xe(p) = {Qe ,X(p) , 
V, := (P^Xeif) = (Qe^Xif), 



■.= P{t,X,t), 
■.= P,{t,X,t), 
V, := iQ,Xeij) = {Pe,X^P) 
:= (P,Xe^) = {Pe,X^P) 
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Vi,x:= J {1- XX)Pdx, Vi,x-= J (1- XX)XPdx, 

Qi,x:= J{l-XX)Qdx, Qi,a:= j{l-XX)XQdx, 

ni^x-= {l-XX)R{x,t)dx, ni^x-= {l-XX)XR{x,t)dx. (5.3.16) 



Note that all quantities in ()5.3.15|) and ()5.3.16p are functions of t alone. Furthermore, 
let 

= e(p{oo) = - / (p{x)dx, = etjj{oo) = - / il){x)dx. (5.3.17) 

Recall from the previous section that when /? = 1 we take to be even. It follows 

that if and ip are even and odd functions respectively. Thus = for (3 = 1, and 
computation gives 

= (7riV)V42-^/-^/^ v^. (5.3.18) 



hm c^ = \l% (5.3.19) 



Hence computation yields 

nS^'oo^^ V 2 
and at t = oo we have 

n,(oo) = 0, g,(oo) = 
Pi,a(oo) = 2c^, Qi,a(oo) = 2c^, 7^i,a(oo) = 0, 
Pi,a(oo) = Qi,a(oo) = :^i,a(oo) = 0. 
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Hence 



'ai,l3i)=Ve, {ai ,(32) = Qe - c^, (ai,/33)=g, 

>3 , /5i 

[a3,P2 
. "3 , /^a 



^ [Pi,A-ai,A(l-^^.)], 

[7^l,A + ai,A (ge - c,^)] 

2 A 

^ [7^i,A + ai,A g^J , 
2 A 



A 

A 
A 

I 
A 



Vi,x - Si.A (1 - v^) 
7^1,A + ai,A - c,^ 
7^i,A + ai,A g, 



(5.3.20) 
(5.3.21) 

(5.3.22) 

(5.3.23) 

(5.3.24) 

(5.3.25) 

(5.3.26) 



As an illustration, let us do the computation that led to ()5.3.22|) in detail. As in 
we use the facts that 5* = S, and {S + S /?*) X = R which can be easily seen by 
writing R = '^'^=i{S X)'' . Furthermore we write R{x,ak) to mean 



lim R(x, y). 

y&J 



In general, since all evaluations are done by taking the limits from within J, we can 
use the identity X 6^ = Sk inside the inner products. Thus 



«2,/92) = Y [{{S + RS) {1 - XX) ,6t - 6^) + ai,x {Qe,St-6^)] 
A 

= t [{{1 - XX) ,{S + R' S) {5t - S^)) + ai,A (Qe(t) - Qe(oo))] 
A 

= t [{{1 - XX) ,{S + R' S)X {6t - 6^)) + ai,A (g. - c, 
A 

= ^ [( (1 - A X) , i?(x, t) - R{x, oo) ) + ai,A (g. - c^)] 
A 

= t [7^l,A(^) - 7^l,A(oo) + ai,A (qe - c^)] 
X 

= Y i'^iA't) + ai,A (ge - c^)] • 

A 
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det {6jk - {aj ,f3k))i 



<j,k<3 



(5.3.27) 



as — s> oo. In order to get our hands on the hmits of the individual terms involved 
in the determinant, we will find differential equations for them first as in [22] • Row 
operation on the matrix show that ai^x and ai^x fall out of the determinant; to see 
this add Aai,A/(2A) times row 1 to row 2 and A Si, a/A times row 1 to row 3. So we 
will not need to find differential equations for them. Our determinant is 



det 



XVi,x 1 A7^l,A A7^lx 



\ 



2A 



1 



A 



2A 

A7^l,A 
A 



1 + 



2A 

A7^l,A 



(5.3.28) 



Proceeding as in j2S] we find the following differential equations 



d 

^^i,A = (A - 7^l,A) , 

^^1,A = Ql,A - A,A, 

^"Pl.A = (^A - 1 - Ui^x 



(5.3.29) 
(5.3.30) 
(5.3.31) 
(5.3.32) 



Let us derive the first equation in (|5.3.3(jp for example. From j22] (equation 2.17), we 
have 

^ = -R{x,t)q^. 
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Therefore 



dt 



d 

It 



Q{x,t)dx-{1- X) [ Q{x,t)dx 



oo 

t 



dQ 



at 



1n+ I -Q^dx-{l-X) 
Qn-In R{x,t)dx - {1 - X)q^ + {1 - X)qj^ R{x,t)dx 

J — oo J oo 

/oo 
(1 - X)R{x,t)dx 
'OO 

^Qn -Qn ^1,A = - ^1,a) • 



Now we change variable from t to s where t = r(s) = 2a\'N+ j^^tb- Then we 
take the hmit — > oo, denoting the hmits of qe,Vi^x, Qi,a, "^i.A? "Pi.A, Qi,a,'^i,a and 
the common hmit of and respectively by q,Vi^x, Qi,x,TZi,x,Vi^x, Qi,x,T^i,x and 
u. We eliminate Qi^x and Qi^x by using the facts that Qi^x = 'Pi,x + X \pl and 
Qi^A = "Pi, A- These limits hold uniformly for bounded s so we can interchange lim 
and Also limAr^oo -^""'^''^Q'jv = ^^"^n-^oo^'^^^Vn — 'i > where q is as in p.l.lO|) . 
We obtain the systems 



d _ 1_ d _ 1 



(5.3.33) 



- ^ (^M - A) , ^7^,, - . (2 P,A + ^2A ) , (5.3.34) 



^nA-^.(l-A-7^,A 



as 



-qV2Vi,x. (5.3.35) 
The change of variables s ^ = J°° q{x, X) dx transforms these systems into constant 
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coefficient ordinary differential equations 



d _ 1 _ d _ 1, _, 
■1-^=^(1^ 1-1 = 7^ , 5.3.36 



. -L in,, _ ,) , jL^., . _L (2 p,, + , (5.3.37) 



Ap„ = (1 - A - n.,) . , = (5.3.38) 

Since lim^^oo = 0, corresponding to the boundary values aXt = oo which we found 
earlier for 'Pi,a, "T^i.a, "Pi.a, 72.i,a, we now have initial values at // = 0. Therefore 

Pi,a(0) = ^i,a(0) = ?i,a(0) = 1i,a(0) = 0. (5.3.39) 



We use this to solve the systems and get 

g = ^^^e'^ + ^^^e-^ (5.3.40) 

u = ^^e^-^^e-^ + K (5.3.41) 

4 4 2' ^ ' 

^i,A = - + A, (5.3.43) 

^i,A = ^ - = (e^ + e-^ - 2). (5.3.44) 

Substituting these expressions into the determinant gives p. 2.11) . namely 



^, ^, ~, A - 1 - cosh/i(s, A) + V A sinh/i(s, A) /r o 

Di(s,A) = D2(s,A) \-2 ^1^^, (5.3.45) 



where Df3 = lim7v-»oo Djs^N- As mentioned in the Introduction, the functional form of 
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the /3 = 1 limiting determinant is very different from what one would expect, unlike 
in the (3 = 4 case. Also noteworthy is the dependence on A = 2A — instead of just 
A. However one should also note that when A is set equal to 1, then A = A = 1. Hence 
in the largest eigenvalue case, where there is no prior differentiation with respect to 
A, and A is just set to 1, a great deal of simplification occurs. The above formula then 
nicely reduces to the /3 = 1 Tracy- Widom distribution. 



Chapter 6 



Applications 



6.1 An Interlacing property 



The following series of lemmas establish Corollary ()1.2.4p : 



Lemma 6.1.1. Define 



A 



Then aj satisfies the following recursion 



A=l 



(6.1.1) 



di = < 



1 ^f J = 0, 

(j - 1) ttj^i for i > 1, j even, 

jttj-i for i >l,i odd. 



(6.1.2) 



Proof. Consider the expansion of the generating function /(A) = \ l around A = 1 



2-A 



j>0 ^' j>0 



Since ttj = j\ bj, the statement of the lemma reduces to proving the following recur- 
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rence for the bi 



Let 



= < 



1 if j = 0, 

^bj_i for j > 1, j even, 
bj-i for j > 1, j odd. 



(6,1,3) 



2- A 



A 



2- A 



2- A 



A 



These are the even and odd parts of / relative to the reflection A — 1 
A — > 2 — A. Recurrence ()6.1.3p is equivalent to 



-(A - 1) or 



which is easily shown to be true. 



□ 



Lemma 6.1.2. Define 



f{s,\) = l 



^ tanh^, 



2- A 



(6.L4) 



forX = 2X- A2. Then 



an 



aA2' 



i2 n+1 



A=l 



2n+ 1 9A2"+i 



A=l 



1 ^/n = 0, 
ifn>l. 



(6.L5) 



Proof. The case n = is readily checked. The main ingredient for the general case is 
Faa di Bruno's formula 



dp" 



9m) = E 
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where k = Y17=i ^« above sum is over all partitions of n, that is all values of 

ki, . . . ,kn such that Yl^=i iK = n. We apply Faa di Bruno's formula to derivatives 



of the function tanh 



J > Ij nonzero only when j = 2, in which case it equals —2. Hence, in 

fl6.1.6p . the only term that survives is the one corresponding to the partition all of 
whose parts equal 2. Thus we have 



2 ' 



which we treat as some function g{X{X)). Notice that for 



n—k 



A=l 



dW^ ^^^^ 



iik = 2j + 1, J > 
for k = 2j, J > 



A=l 



tanh ^^ ^ 



Q ^2n+l-fc 



A=l 



tanh 



iik = 2j, J > 

for k = 2j + 1, J >0 



A=l 



Therefore, recalling the definition of aj in (jG.l.lj) and setting k = 2j, we obtain 



i2n 



2n 



f{s,X) 



A=l 



2n 

E 

fc=0 



2n\ (9^ 



k dX'' y 2-XdX 



2n-k 



A=l 



^ (2n)!(-l)"-^' 9"-^' , /i(s,A) 



j=0 



{2j)\{n-j)\ dX^~i 



A=l 
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Similarly, using k = 2 j + 1 instead yields 



g2 n+l 



/(^,A) 



A=l 



2 n+l 



k=0 



(2 72 + 1) ^ 



i=o 

^2„ 



(2j)!(^-j)! 2j + l dX'^-i 



tanh 



A=l 

/i(s, A) 



A=l 



(2^ + l)^/(^'^) 



A=l 



since .,^/(2 j + 1) = . Rearranging this last equality leads to ()6.1.5p . 



□ 



Lemma 6.1.3. LetDi{s,X) and Di{s,~X) he as m ffTTT]) and (frT2|) . T/ien 



D^{s,\) = D^{s,\) 1 



^ tanh/^^^ 



2- A 



(6.1.7) 



Proof. Using the facts that —1 — coshx = —2 cosh^ |, 1 = cosh^x — sinh^x and 



sinh X = 2 sinh | cosh | we get 



-1 -cosh/i(g,A) u^n^ ~ A+VAsinh/i(s,A) 

-Di(s, A) = r D-2[s,\) + D2[s,\) — 



A-2 



A-2 



A-2 



+D2{s, A) 



A (cosh^ ^ - sinh^ ^) + sinh/i(s, A) 



A-2 



A(s,A) 



^4(5, A) A sinh^ - \/I sinh /i(s, A) 



cosh 



2 / 



2- A 



' X sinh' - 2 sinh ( cosh f ^) 
(2 - A) cosh' (h(£^ 



Z^4(^, A) I 1 - 2 ,/ — tanh + -A_ tanh' 



2- A 



A(^,A) tanh^^) 



2 2 - A 

2 
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□ 

For notational convenience, define di{s, A) = /^^^(s, A), and ^4(5, A) = /^^^(s, A). 
Then 

Lemma 6.1.4. For n > 0, 

1 Q2n+1 gin 

_~(2r2 + 1)! 9A2"+i ^ {2n)\ dl^ 



di{s, A) 



-1)" 9" 



A=l 



9 A" 



6/4(5, A) 



Proof. Let 



/(5,A) = 1 



^ tanh^ 



2- A 



by the previous lemma, we need to show that 



g2 n+l 



1 d' 



2n 



(2n+ 1)! 9A2"+i (2n)! SA^' 



^(s. A) /(s. A) 



A=l 



(6.L 



equals 



^ ^^4(5, A) 



9A 

Now formula ()6.1.6p applied to ^4(5, A) gives 



A=l 



Qk 



: d4^{s, A) 



A=l 



if fc = 2j + l, j>0, 

^^4(s,A) iffc = 2j,j>0. 
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Statistic 


Eigenvalue 




a 


7i 


72 


Ai 


-1.206548 


1.267941 


0.293115 


0.163186 


A2 


-3.262424 


1.017574 


0.165531 


0.049262 


A3 


-4.821636 


0.906849 


0.117557 


0.019506 


A4 


-6.162036 


0.838537 


0.092305 


0.007802 



Table 6.1: Mean, standard deviation, skewness and kurtosis data for first four edge- 
scaled eigenvalues in the 13 = 1 Gaussian ensemble. Corollarv II .2.41 implies that rows 
2 and 4 give corresponding data for the two largest eigenvalues in the (3 = 4 Gaussian 
ensemble. Compare to Table 1 in [53], keeping in mind that the discrepancy in the 
(3 = 4 data is caused by the different normalization in the definition of -^4(5, 1); see 
footnote in the Introduction, as well as the comments following Equation ()4.3.4H1 . 
and Section IX!4l 



Therefore 



2 n+l 



(2n + l)! aA2'^+i 



c?4(s, A)/(s, A) 



A=l 
2 n+l 



(2n + 

n 



1)! ^ 



fc=0 

:-iv 



2n+l\ , Q2n+l^k 

a. 



A=l 



Qj Q2n-2j+l 



3=0 



A=l 



Similarly, 



d' 



,2n 



(2n)! dX 



2n 



di{s,X) f{s,X) 



E 

3=0 



{2n-2j)\j\ dXi ^A^'^-^J 



A=l 



A=l 



Therefore the expression in ()6.1.8p equals 



3=0 



.1)3 



(2n-2j)!j! dXi 



— (s. A) 



Q2n-2j 



f- 



1 



^2n-2j+l 



dX2n~2j-' 2n-2j + 1 (9A2"-2j+i 



/ 



A=l 
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Statistic 


Eigenvalue 




a 


7i 


72 


Ai 


-1.771087 


0.901773 


0.224084 


0.093448 


A2 


-3.675440 


0.735214 


0.125000 


0.021650 



Table 6.2: Mean, standard deviation, skewness and kurtosis data for first two edge- 
scaled eigenvalues in the (3 = 2 Gaussian ensemble. Compare to Table 1 in |25j . 



Now Lemma f6 . 1 . 21 shows that the square bracket inside the summation is zero unless 
j = n, in which case it is 1. The result follows. □ 

In an inductive proof of Corollarv ll.2.4[ the base case ^4(5, 2) = 1) is easily 
checked by direct calculation. Lemma [6 . 1 .41 establishes the inductive step in the proof 
since, with the assumption F4(s,n) = Fi{s, 2n), it is equivalent to the statement 

F4(s,n + 1) = Fi{s,2n + 2). 



6.2 The Wishart ensembles 

Consider an n x p data matrix X whose rows Xi are independent Gaussian A'j,(0, S). 
The product X, of X and its transpose is (up to factor 1/n) a sample covariance 
matrix. In this setting, the matrices A = X^ X are said to have Wishart distribution 
Wp{n, S). The so-called "Null Case" corresponds to S = Id„. Let Ai > ■ • • > Ap be 
eigenvalues of X* X and define 

Johnstone proved the following theorem (see [TH]). 
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Theorem 6.2.1. Ifn/p — > 7 > then 

^1 - f^np 



a 

np 



Fi{s,l). 



Soshnikov generalized this to the m*^ largest eigenvalue (see|19j). 
Theorem 6.2.2. Ifn/p ^ 7 > then 



a 

np 



Fi{s,m). 



Note that Soshnikov redefines fi^^, a^^ by letting n ^ n + lin Johnstone definition, 
but this does not affect the limiting distribution. El Karoui proves lU'.2.1l for < 7 < 00 
(see jin])- It is an open problem to show Soshnikov's theorem 16.2.21 holds in this 
generality, but if it does then our results will also have wider applicability. Soshnikov 
went further and disposed of the Gaussian assumption, showing the universality of 
this result. Let us redefine the matrices X = {xij} to satisfy 

• Exij = 0, E{xijY = 1- 

• the r.v's Xij have symmetric laws of distribution 

• all moments of these r.v's are finite 

• the distributions of the Xij decay at least as fast as a Gaussian at infinity: 
E(xij)2"^ < (const m)"' 

. n-p = 0(pi/3) 
With these new assumptions, Soshnikov proved the following fsee|19j). 
Theorem 6.2.3. Ifn/p — * 7 > then 



a 

np 



Fi(s,m). 
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100 X 100 


100 X 400 


i^i-Percentile 


M 




^3 


M 




^3 


U.Ui 


U.UUo 


u.uuo 


n nn/1 

U.UU4I: 


n nns 

U.UUo 


U.UUO 


n nn/1 

U.UU4I: 


U.Uo 


U.U4Z 


U.Uoo 


n not; 
u.uzo 


U.U4Z 


U.Uo ( 


n n'?o 


U.iU 


n nan 
u.uyu 


U.U ( o 


u.uoy 


U.Uoo 


n nsi 

U.Uoi 


n nfiR 

U.UOO 


U.oU 


n 9Q/1 


U.ZOo 


U.ZOO 


n OS'? 

U.ZOO 


n 9R7 


n 9c;/i 

U.Z04I: 


U.OU 


n zLQ? 


n 477 


n zLzLD 




n 471 


n 


U. ( u 


0.699 


0.690 


0.659 


0.685 


0.679 


0.669 


0.90 


0.902 


0.891 


0.901 


0.898 


0.894 


0.884 


0.95 


0.951 


0.948 


0.950 


0.947 


0.950 


0.941 


0.99 


0.992 


0.991 


0.991 


0.989 


0.991 


0.989 



Table 6.3: Percentile comparison of Fi vs. empirical distributions for 100 x 100 and 
100 X 400 Wishart matrices with identity covariance. 

Our results are immediately applicable here and thus explicitly give the distribu- 
tion of the m^^ largest eigenvalue of the appropriate Wishart distribution. Table 16.31 
shows a comparison of percentiles of the Fi distribution to corresponding percentiles 
of empirical Wishart distributions. Here Aj denotes the i^^ largest eigenvalue in the 
Wishart Ensemble. The percentiles in the Aj columns were obtained by finding the 
ordinates corresponding to the Fi-percentiles listed in the first column, and comput- 
ing the proportion of eigenvalues lying to the left of that ordinate in the empirical 
distributions for the Aj. The bold entries correspond to the levels of confidence most 
commonly used in statistical applications. The reader should compare this table to 
a similar one in [T5] . 



Chapter 7 



Numerics 



7.1 Partial derivatives of q{x,X) 

Let 

Qn{x) = -^q{x,X) , (7.1.1) 

so that go equals q from p.l.8|) . In order to compute F/^{s, m) it is crucial to know qn 
with < 77, < m accurately. Asymptotic expansions for g„ at — oo are given in j22j . 
We outline how to compute go and gi as an illustration. From [22], we know that, as 
t +00, go(— 1/2) is given by 

1 ^ / 1 73 10657 13912277 I 
-v^ 1 - - - — - — ^ + O — 



t3 2t6 2^9 8ti2 Vt^s 

whereas gi(— 1/2) can be expanded as 



exp(k3/2) / 17 1513 850193 407117521 ^ 

1 + TTT^ + TTTTTTr^ + . ^.r^r.. +0 



2v^tV4 V 24t3/2 ^ 2^32^3 ^ 21034^9/2 21535^6 ^ U15/2 



87 
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Quantities needed to compute Ffs{s,m),m = 1,2, are not only go and qi but also 
integrals involving qq, such as 



(x — s) go(x) dx, Jo 



qo{x) dx. 



(7.1.2) 



Instead of computing these integrals afterward, it is better to include them as variables 
in a system together with go, as suggested in [TH]. Therefore all quantities needed 
are computed in one step, greatly reducing errors, and taking full advantage of the 
powerful numerical tools in MATLAB"'"'^ . Since 



ql{x)dx, Io=ql, Jq = -Qo, 



(7.1.3) 



the system closes, and can be concisely written 



d 

ds 



\Jo J 



sqo + 2ql 



\ 



^1 



-qo 



J 



(7.1.4) 



We first use the MATLAB™ built-in Runge-Kutta-based ODE solver ode45 to 
obtain a first approximation to the solution of ()7.1.4p between x = 6, and x = —8, 
with an initial values obtained using the Airy function on the right hand side. Note 
that it is not possible to extend the range to the left due to the high instability 
of the solution a little after —8; (This is where the transition region between the 
three different regimes in the so-called "connection problem" lies. We circumvent 
this limitation by patching up our solution with the asymptotic expansion to the left 
of X = —8.). The approximation obtained is then used as a trial solution in the 
MATLAB™ boundary value problem solver bvp4c, resulting in an accurate solution 
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vector between x = 6 and x = —10. 
Similarly, if we define 



{x — s) qo{x) qi{x) dx, Ji 



then we have the first-order system 



d 

ds 



h 

V ^1 J 



( 



qoix) qi{x)dx, 



\ 



sqi + Qql qi 

Qoqi 
-qoQi 



J 



(7.1.5) 



(7.1.6) 



which can be implemented using bvp4c together with a "seed" solution obtained in 
the same way as for go- 



7.2 MATLAB^^^^ code 

In this section we gather the MATLAB™ code used to evaluate and plot the distri- 
butions. It is broken down into twenty- five functions/files for easier use. They are 
implemented in version 7 of MATLAB™ , and may not be compatible with earlier 
versions. The function "twmakenew" should be ran first. It initializes everything 
and creates some data files which will be needed. There will be a few warning mes- 
sages that can be safely ignored. Once that is done, the distributions are ready to be 
used. The only functions directly needed to evaluate them are "twdens" and "twdist" 
which, as their names indicate, give the Tracy- Widom density and distribution func- 
tions. They both take 3 arguments; the first is "beta" which is the beta of RMT so 
it can be 1, 2 or 4; then "n" which is the eigenvalue need; finally "s" which is the 
value where you want to evaluate the function. All these functions are vectorized. 
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meaning that they can take a vector as well as scalar argument. This is convenient for 
plotting. One does not need to run "twmakenew" every time you want to work with 
these functions. In fact one should only need to run it once ever, unless some files are 
deleted or moved. Once it is ran the first time, it creates data files that contain all 
you need to evaluate the functions. On subsequent uses one can just start by running 
"tw" instead, which will initialize by reading all the data that "twmakenew" produced 
the first time it ran. This is much faster. The only reason for running "twmakenew" 
again is to clear everything and refresh all the data. Otherwise one can just start 
with "tw" and then start using "twdens" and "twdist" . As a summary here is a list 
of commands that will plot the density of the S^'^ GOE eigenvalue between and —9: 

twmakenew 
s=-9:0.1:0 

plot(s, twdens (1, 3, s)) 

The MATLAB^^ functions needed are below. Note that the digits at the begin- 
ning of each line are not part of the code and are added only for (ease of) reference 
purposes. They should be omitted when coding. These functions are available elec- 
tronically by writing to the author. 



CHAPTER 7. NUMERICS 91 
7. 

7. (c) 2004-2005 Momar Dieng, All rights reserved. 

7. 

y„ Comments? e-mail momar@math.ucdavis.edu 

°/o This software was developed with support of the 

7. National Science Foundation/ Grant DMS— 0304414 



In EVERY case, the software is COPYRIGHTED by the original author. 
7o For permissions, see below. 

7. 

y„ Comments/Permissions/Information: e-mail 

7o momarOmath . ucdavis . edu or tracySmath . ucdavis . edu 

y. 

y. COPYLEFT NOTICE. 

% Permission is granted to make and distribute verbatim copies of this 
% entire software package, provided that all files are copied 
y **together** as a **unit** . 

y. 

y Here **copying as a unit** means: all the files listed are copied, 
y Permission is granted to make and distribute modified copies of this 
y software, under the conditions for verbatim copying, provided that 
% the entire resulting derived work is distributed under the terms 
% of a permission notice identical to this one. Names of new authors, 
y their affiliations and information about their improvements may be 
% added to the files. 

y 

y The purpose of this permission notice is to allow you to make 
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% copies of the software and distribute them to others, for free 

% or for a fee, subject to the constraint that you maintain the 

% collection of tools here as a unit. This enables people you 

% give the software to to be aware of its origins, to ask questions 

y„ of us by e-mail, to request improvements, obtain later releases, 

°/„ etc . . . 

% 

% If you seek permission to copy and distribute translations of 
y„ this software into another language, please e-mail a specific 
% request to momar@math.ucdavis.edu or tracy@math.ucdavis.edu. 
% 

% If you seek permission to excerpt a **part** of the software for 

% example to appear in a scientific publication, please e-mail a 

% specific request to momar@math.ucdavis.edu or tracy@math.ucdavis.edu. 

y. 



1 function [dbcdya,dbcdyb] = qObcjac(ya,yb) 



2 


dbcdya = [ 


1 














3 


















4 


















5 


















6 















0] 


7 


dbcdyb = [ 

















8 




1 














9 










1 








10 













1 





11 
















1]; 
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1 function res = qObc(ya,yb) 

2 global qOsolbcleft; 

3 global qOsolbcright ; 

4 res = [ya(l) - qOlef texp(-qOsolbclef t) ; yb(l) - . 

5 airy (qOsolbcright) ; yb(3); yb(4); yb(5)] ; 

1 function y = qObvpinit(s) 

2 y = [ qO(s) 

3 qOprime(s) 

4 xmuO(l,s)-s*muO(l,s) 

5 -mu0(2,s) 

6 muO(l,s) 

7 ]; 

1 function jac = qOjac(x,y,e) 

2 jac = [ 10 

3 x+6*y(l)-2 

4 1 

5 2*y(l) 

6 -1 0]; 

1 function y = qOleftexp(x) 

2 y = sqrt(x./2) .*(1 - l./(8*x."3) - 73 . / (128*x . ~6) 

3 10657./(1024*x."9) - 13912277. /(32768*x. "12) 

4 - (8045883943/262144) ./x." 15 - ... 

5 (14518451390349/4194304) ./x.~ 18 - ... 



CHAPTER 7. NUMERICS 



6 (18847128706420641/33554432) ./x. "21 - ... 

7 (266287398541797779277/2147483648) ./x.~24 ... 

8 - (614077537500104697967243/17179869184) ./x.~27 

9 - (3583820408994355704357070999/ . . . 

10 274877906944) ./x. "30 - ... 

11 (12909064475999593160583566233991/ . . . 

12 2199023255552) ./x. "33 - ... 

13 (225032454118493801360210306128796801/ . . . 

14 70368744177664) ./x. "36 - ... 

15 (1167024146596038918316117313883961507459/ . . . 

16 562949953421312) ./x. "39 - ... 

17 (14205068353420952020935752835796596130983637/ 

18 9007199254740992) ./x. "42) ; 



1 function y = qOleftexp(x) 

2 y = l./4.*2."(l./2) ./x."(l./2) .*(l-l./8./x."3-73./128./. . . 

3 x. "6-10657. /1024./X. "9-13912277. /32768./. . . 

4 X. "12-8045883943. /262144./X. "15-14518451390349./. . . 

5 4194304. /x." 18-588972772075645. /1048576./. . . 

6 X . "21-8126446488702325 . /65536 . /x . "24-4575234185905045 . / . . . 

7 128 . /x . "27-13037862696344220 . /x . "30-5870362872883376128 . / . . . 

8 X. "33-3197903511683275948032. /x. "36-. . . 

9 2073051324551114335977472. /x. "39-. . . 

10 1577079395234210168300371968. /x. "42) +1. /2. *2. " (1. /. . . 

11 2) .*x."(l./2) .* (3. /8./x. "4+219. /64./X. "7+95913. /1024./. . . 

12 x."10+41736831./8192./x. "13+120688259145. /262144./X. "16+. . . 

13 130666062513141 . /2097152 . /x . " 19+12368428213588545 . /1048576 . / . 
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14 



X . "22+24379339466106975 . /8192 . /x . "25+123531323019436215 . / . . . 



15 



128 . /x . "28+391135880890326600 . /x . "31+193721974805151412224 . / . 



16 



X. "34+115124526420597934129152. /x. "37+. . . 



17 



80849001657493459103121408. /x. "40+. . . 



18 



66237334599836827068615622656 . /x . "43) ; 



1 function y = qO(x) 

2 global qOsol; 

3 global qOsolbcleft; 

4 global qOsolbcright ; 

5 y = zeros(size(x) ) ; 

6 ix = x < qOsolbcleft; 

7 if (sumdx)) 

8 y(ix) = q01eftexp(-x(ix) ) ; 

9 end 

10 ix = qOsolbcleft <= x & x <= qOsolbcright; 

11 if (suni(ix)) 

12 y(ix) = deval (qOsol ,x(ix) , 1) ; 

13 end 

14 ix = x > qOsolbcright; 

15 if (suin(ix)) 

16 y(ix) = airy(x(ix)); 

17 end 



1 function dydx = qOode(x,y) 

2 dydx = [y(2); x*y(l)+2*y(l) "3; y(4); y(l)"2; -y(l)] ; 
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1 function y = qOprime(x) 

2 prec = 1/100; 

3 if (max(x)-min(x) > 10*prec) 

4 linspace = min(x) :prec :max(x) ; 

5 qOspline = spline (linspace, qO (linspace) ) ; 

6 y = fnval (fnder (qOspline) , x) ; 

7 else 

8 linspace = min(x)-l :prec :max(x)+l ; 

9 qOspline = spline (linspace, qO (linspace) ) ; 

10 y = fnval (fnder (qOspline) , x) ; 

11 end 

1 function y = qOseed(s) 

2 global qOguessbclef t ; 

3 global qOguessbcright ; 

4 global qOguess; 

5 is = s < qOguessbclef t ; 

6 if (sum(is)) 

7 y(:,is) = [yl (s (is) , qOguessbclef t , qOguessbcright) ; . 

8 y2 (s (is) , qOguessbclef t, qOguessbcright) ; . 

9 y3(s(is), qOguessbclef t, qOguessbcright) ; 

10 y4 (s (is) , qOguessbclef t, qOguessbcright) ; . 

11 y5 (s (is) , qOguessbclef t, qOguessbcright) ] 

12 end 

13 is = qOguessbclef t <= s & s <= qOguessbcright; 

14 if (suni(is)) 

15 y(:,is) = deval (qOguess , s (is) ) ; 
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16 end 

17 is = s > qOguessbcright ; 

18 if (suni(is)) 

19 y(:,is) = [airy (s (is) ) ; airy (1 , s (is) ) ; ; airy (s (is) ) "2 ; 0] ; 

20 end 

21 function y = yl(x,a,b) 

22 y = qOleftexp(-x) ; 

23 function y = y2(x,a,b) 

24 y = qOleftexpprime(-x) ; ; 

25 function y = y3(x,a,b) 

26 global qOguess; 

27 y = zeros (size (x) ) ; 

28 for i=l:length(x) 

29 y(i) = quad('x.*q01eftexp(-x) ."2\x(i) ,a,le-14) + ... 

30 deval(q0guess,a,3) - a*deval(qOguess,a,4) + . . . 

31 x(i)*y4(x(i),a,b); 

32 end 

33 function y = y4(x,a,b) 

34 global qOguess; 

35 y = zeros (size (x) ) ; 

36 for i=l:length(x) 

37 y(i) = quad('-q01eftexp(-x) ."2^x(i) ,a,le-14) + ... 

38 deval(q0guess,a,4) ; 

39 end 

40 function y = y5(x,a,b) 

41 global qOguess; 

42 y = zeros (size (x) ) ; 
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43 for i=l:length(x) 

44 y(i) = quad('q01eftexp(-x) ^x(i) ,a,le-14) + 

45 deval (qOguess , a, 5) ; 

46 end 



1 function [dbcdya,dbcdyb] = qlbcjac (ya,yb) 



2 


dbcdya = [ 
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dbcdyb = [ 
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1] 





1 function res = qlbc(ya,yb) 

2 global qlsolbcleft; 

3 global qlsolbcright ; 

4 res = [ya(l) - qllef texp(-qlsolbclef t) ; yb(l) - (1/2)* 

5 airy (qlsolbcright) ; yb(3); yb(4); yb(5)] ; 



1 function y = qlbvpinit(s) 

2 y = [ ql(s) 

3 qlprime(s) 
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4 xmuOmul(s)-s*muOmul(s) 

5 -muOmul(s) 

6 mul(l,s) 

7 ]; 

1 function jac = qljac(x,y) 

2 jac = [ 10 

3 x+6*q0(x)~2 

4 1 

5 q0(x) 

6 -1 0]; 

1 function y = qlleftexp(x) 

2 y= (l/(2*sqrt(2*pi)))*(exp(((2*x) .-(3/2))/3) ./. 

3 ((2*x) .-(1/4))) .*(1 + 17./(24*(2*x) .-(3/2)) 

4 + 1513./((2-7)*(3~2)*(2*x) .~3) + 850193./.. 

5 ((2~10)*(3-4)*(2*x) .-(9/2)) +407117521./.. 

6 ((2~15)*(3-5)*(2*x) .~6)) ; 



1 function y = ql(x) 

2 global qlsol; 

3 global qlsolbcleft; 

4 global qlsolbcright ; 

5 y = zeros(size(x) ) ; 

6 ix = X < qlsolbcleft; 

7 if (sum(ix)) 
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8 y(ix) = qllef texp(-x(ix) ) ; 

9 end 

10 ix = qlsolbcleft <= x & x <= qlsolbcright ; 

11 if (suni(ix)) 

12 y(ix) = deval(qlsol,x(ix) , 1) ; 

13 end 

14 ix = X > qlsolbcright; 

15 if (sum(ix)) 

16 y(ix) = (l/2)*airy(x(ix)) ; 

17 end 

1 function dydx = qlode(x,y) 

2 dydx = [y(2); x*y(l) + 6* (qO (x) "2) *y (1) ; y(4) 

3 qO(x)*y(l); -y(l)] ; 

1 function y = qlprime(x) 

2 prec = 1/100; 

3 if (max(x)-min(x) > 10*prec) 

4 linspace = min(x) :prec :max(x) ; 

5 qlspline = spline (linspace, ql (linspace) ) ; 

6 y = fnvaKfnder (qlspline) , x) ; 

7 else 

8 linspace = min(x)-l :prec :max(x)+l ; 

9 qlspline = spline (linspace, ql (linspace) ) ; 

10 y = fnval (fnder (qlspline) , x) ; 

11 end 
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1 function y = qlseed(x) 

2 global qlguess; 

3 y = deval (qlguess ,x) ; 



1 % The function twdens(beta, n, s) computes the 

2 y„ (edge-scaled) limiting density for the n"th 

3 7o largest eigenvalue 

4 y„ of the specified beta-ensemble at s. 

5 function y = twdens(beta, n, s) 

6 global qOsol; 

7 global qOsolbcleft; 

8 global qOsolbcright ; 

9 global qlsol; 

10 global qlsolbcleft; 

11 global qlsolbcright ; 

12 if (beta == 1) 

13 if (n == 1) 

14 is = s < qOsolbcleft; 

15 if (sum(is)) 

16 yds) = 0; 

17 end 

18 is = qOsolbcleft <= s & s <= qOsolbcright; 

19 if (sum(is)) 

20 y(is) = (-deval(qOsol,s(is) ,4)+ ... 

21 qO(s(is))) .*twdist(l, l,s(is))/2; 

22 end 

23 is = s > qOsolbcright; 
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24 if (sum(is)) 

25 y(is) = 0; 

26 end 

27 elseif (n == 2) 

28 is = s < qOsolbcleft; 

29 if (sum(is)) 

30 y(is) = 0; 

31 end 

32 is = qOsolbcleft <= s & s <= qOsolbcright ; 

33 if (sum(is)) 

34 y(is) = (-deval(qOsol,s(is) ,4) - tanh(deval( . . . 

35 qOsol,s(is) ,5)/2) .*qO(s(is))) . *twdist (1 , 2 , s (is) ) /2 ; 

36 end 

37 is = s > qOsolbcright; 

38 if (sum(is)) 

39 y(is) = 0; 

40 end 

41 elseif (n == 3) 

42 is = s < qOsolbcleft; 

43 if (sum(is)) 

44 y(is) = 0; 

45 end 

46 is = qOsolbcleft <= s & s <= qOsolbcright; 

47 if (sum(is)) 

48 y(is) = twdensd, l,s(is)) .*(ones(size(s(is))) ... 

49 - exp(deval(qOsol,s(is) ,5)) + ... 

50 4*deval(qlsol,s(is) ,3) + ... 
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51 2*deval(qlsol,s(is) ,5))/4 + ... 

52 twdist(l,l,s(is)) .*(exp(deval( . . 

53 qOsol,s(is) ,5)) .*qO(s(is)) + ... 

54 4*deval(qlsol,s(is) ,4) - .. 

55 2*ql(s(is)))/4 + twdens (1 , 2 , s (is) ) ; 

56 end 

57 is = s > qOsolbcright ; 

58 if (sum(is)) 

59 y(is) = 0; 

60 end 

61 elseif (n == 4) 

62 is = s < qOsolbcleft; 

63 if (sum(is)) 

64 y(is) = 0; 

65 end 

66 is = qOsolbcleft <= s & s <= qOsolbcright; 

67 if (sum(is)) 

68 y(is) = twdens (1 , 1 , s(is) ). *exp(deval (.. . 

69 qOsol,s(is) ,5)/2) .*((ones(size(s(is))) ... 

70 + 2*deval(qlsol,s(is) ,3)) .*sinh(deval(. . . 

71 q0sol,s(is) ,5)/2) - cosh(deval( . . . 

72 q0sol,s(is) ,5)/2) .*deval(qlsol,s(is) ,5))/. . . 

73 2 + (twdist(l,l,s(is)) .*exp(deval(. . . 

74 q0sol,s(is) ,5) /2) . * (-qO(s (is) )/2) . * ( (ones ( . . . 

75 size(s(is))) + 2*deval(qlsol,s(is) ,3)) .*sinh( . 

76 deval(qOsol,s(is) ,5)/2) - cosh(deval( . . . 

77 qOsol,s(is) ,5)/2) .*deval(qlsol,s(is) ,5)))/2 + 
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78 (twdist(l,l,s(is)) . *exp(deval(qOsol , s(is) ,5)/. . . 

79 2) .*(2*deval(qlsol,s(is) ,4) . *sinh(deval( . . . 

80 qOsol,s(is) ,5)/2) + (ones(size(s(is))) + ... 

81 2*deval(qlsol,s(is) ,3)) . *cosh(deval( . . . 

82 qOsol,s(is) ,5)/2) .*(-q0(s(is))/2) - ... 

83 sinh(deval(qOsol,s(is) ,5)/2) .*(-qO(s(. . . 

84 is))/2) .*deval(qlsol,s(is) ,5) +cosh(... 

85 deval(qOsol,s(is) , 5)/2) .*ql(s(is))))/2 + ... 

86 twdens(l,3,s(is)) ; 

87 end 

88 is = s > qOsolbcright ; 

89 if (sum(is)) 

90 y(is) = 0; 

91 end 

92 else 

93 fprintf ( ['Don' 't know how to plot the density for ' ... 

94 ['that GOE eigenvalue. Check that the ' ... 

95 'data for it is available.']); 

96 end 

97 elseif (beta == 2) 

98 if (n == 1) 

99 is = s < qOsolbcleft; 

100 if (sum(is)) 

101 y(is) = 0; 

102 end 

103 is = qOsolbcleft <= s & s <= qOsolbcright; 

104 if (sum(is)) 
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105 y(is) = twdist(2, l,s(is)) .*(-deval(qOsol,s(is) ,4)) ; 

106 end 

107 is = s > qOsolbcright ; 

108 if (sum(is)) 

109 y(is) = 0; 

110 end 

111 elseifCn == 2) 

112 is = s < qOsolbcleft; 

113 if (sum(is)) 

114 y(is) = 0; 

115 end 

116 is = qOsolbcleft <= s & s <= qOsolbcright; 

117 if (sum(is)) 

118 y(is) = twdens(2, 1, s(is)) .*(2*deval(qlsol,s(is) ,3) 

119 + ones(size(s(is)))) - twdist(... 

120 2, 1, s(is)) .*(2*(-deval(qlsol,s(is) ,4))) ; 

121 end 

122 is = s > qOsolbcright; 

123 if (sum(is)) 

124 y(is) = 0; 

125 end 

126 else 

127 fprintf ( ['Don' 't know how to plot the density for ' ... 

128 ['that CUE eigenvalue. Check that the data ' . 

129 'for it is available.']); 

130 end 



131 elseif (beta == 4) 
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132 if (n == 1) 

133 is = s < qOsolbcleft; 

134 if (sum(is)) 

135 y(is) = 0; 

136 end 

137 is = qOsolbcleft <= s & s <= qOsolbcright ; 

138 if (sum(is)) 

139 y(is) = (l./(2*sqrt(twdist(2,l,s(is))))) .*twdens(. 

140 2,l,s(is)) .*cosh(deval(qOsol,s(is) ,5)/2) + ... 

141 sqrt(twdist(2,l,s(is))) . *sinh(deval ( . . . 

142 qOsol,s(is) ,5)/2) . * (-qO (s (is) ) /2) ; 

143 end 

144 is = s > qOsolbcright; 

145 if (sum(is)) 

146 y(is) = 0; 

147 end 

148 elseif (n == 2) 

149 is = s < qOsolbcleft; 

150 if (sum(is)) 

151 y(is) = 0; 

152 end 

153 is = qOsolbcleft <= s & s <= qOsolbcright; 

154 if (sum(is)) 

155 y(is) = twdens(l,4,s(is)) ; 

156 end 

157 is = s > qOsolbcright; 

158 if (sum(is)) 
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159 y(is) = 0; 

160 end 

161 else 

162 fprintf ( ['Don' 't know how to plot the density for that ' 

163 ['GSE eigenvalue. Check that the data for it ' 

164 'is available. ']) ; 

165 end 

166 else 

167 fprintf ( 'Don' 't know anything about that beta!'); 

168 end 



1 % The function twdist(beta, n, s) computes the 

2 % (edge-scaled) limiting distribution for the 

3 % n"th largest eigenvalue 

4 % of the specified beta-ensemble at s. 

5 function y = twdist(beta, n, s) 

6 tw; 

7 global qOsol; 

8 global qOsolbcleft; 

9 global qOsolbcright ; 

10 global qlsol; 

11 global qlsolbcleft; 

12 global qlsolbcright ; 

13 if (beta == 1) 

14 if (n == 1) 

15 is = s < qOsolbcleft; 

16 if (sum(is)) 
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17 y(is) = 0; 

18 end 

19 is = qOsolbcleft <= s & s <= qOsolbcright ; 

20 if (sum(is)) 

21 y(is) = sqrt(twdist(2,l,s(is))) .*exp(. . . 

22 -deval (qOsol , s (is) , 5) /2) ; 

23 end 

24 is = s > qOsolbcright; 

25 if (sum(is)) 

26 y(is) = 1; 

27 end 

28 elseif (n == 2) 

29 is = s < qOsolbcleft; 

30 if (sum(is)) 

31 y(is) = 0; 

32 end 

33 is = qOsolbcleft <= s & s <= qOsolbcright; 

34 if (sum(is)) 

35 y(is) = sqrt(twdist(2, l,s(is))) .*sinh(deval( . 

36 qOsol,s(is) ,5)/2) + twdist (1 , 1 , s (is) ) ; 

37 end 

38 is = s > qOsolbcright; 

39 if (sum(is)) 

40 yds) = 1; 

41 end 

42 elseif (n == 3) 

43 is = s < qOsolbcleft; 
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44 if (sum(is)) 

45 y(is) = 0; 

46 end 

47 is = qOsolbcleft <= s & s <= qOsolbcright ; 

48 if (sum(is)) 

49 y(is) = (l/4)*twdist(l,l,s(is)) .*(ones(size(s(is))) . 

50 - exp(deval(qOsol,s(is) ,5)) + 4*deval( . . . 

51 qlsol,s(is) ,3) + 2*deval(qlsol,s( . . . 

52 is), 5)) + twdist(l,2,s(is)) ; 

53 end 

54 is = s > qOsolbcright; 

55 if (sum(is)) 

56 y(is) = 1; 

57 end 

58 elseif (n == 4) 

59 is = s < qOsolbcleft; 

60 if (sum(is)) 

61 y(is) = 0; 

62 end 

63 is = qOsolbcleft <= s & s <= qOsolbcright; 

64 if (sum(is)) 

65 y(is) = (l/2)*twdist(l,l,s(is)) .*exp(deval(qOsol,s(. 

66 is) ,5)/2) .*((ones(size(s(is))) + 2*deval(qlsol, s(. 

67 is) ,3)) .*sinh(deval(qOsol,s(is) ,5)/2) - cosh(deval( 

68 qOsol,s(is) ,5)/2) .*deval(qlsol,s(is) ,5)) + ... 

69 twdist(l,3,s(is)) ; 

70 end 
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71 is = s > qOsolbcright ; 

72 if (sum(is)) 

73 y(is) = 1; 

74 end 

75 else 

76 f printf ( [ 'Don' know how to plot the distribution 

77 ['for that GOE eigenvalue. Check that the 

78 'data for it is available.']); 

79 end 

80 elseif (beta == 2) 

81 if (n == 1) 

82 is = s < qOsolbcleft; 

83 if (sum(is)) 

84 y(is) = 0; 

85 end 

86 is = qOsolbcleft <= s & s <= qOsolbcright; 

87 if (sum(is)) 

88 y(is) = exp(-deval (qOsol , s(is),3)); 

89 end 

90 is = s > qOsolbcright; 

91 if (sum(is)) 

92 y(is) = 1; 

93 end 

94 elseif (n == 2) 

95 is = s < qOsolbcleft; 

96 if (sum(is)) 

97 y(is) = 0; 
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98 end 

99 is = qOsolbcleft <= s & s <= qOsolbcright ; 

100 if (sum(is)) 

101 y(is) = twdist(2,l,s(is)) .*(2*deval(qlsol, s(is),3)... 

102 + ones(size(s(is)))) ; 

103 end 

104 is = s > qOsolbcright; 

105 if (sum(is)) 

106 y(is) = 1; 

107 end 

108 else 

109 fprintf ( ['Don' know how to plot the distribution ' ... 

110 ['for that GUE eigenvalue. Check that the ' ... 

111 'data for it is available.']); 

112 end 

113 elseif (beta == 4) 

114 if (n == 1) 

115 is = s < qOsolbcleft; 

116 if (sum(is)) 

117 y(is) = 0; 

118 end 

119 is = qOsolbcleft <= s & s <= qOsolbcright; 

120 if (sum(is)) 

121 y(is) = sqrt(twdist(2, l,s(is))) .*cosh(deval( . . . 

122 qOsol,s(is) ,5)/2) ; 

123 end 

124 is = s > qOsolbcright; 
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125 if (sum(is)) 

126 y(is) = 1; 

127 end 

128 elseif (n == 2) 

129 is = s < qOsolbcleft; 

130 if (sum(is)) 

131 y(is) = 0; 

132 end 

133 is = qOsolbcleft <= s & s <= qOsolbcright ; 

134 if (sum(is)) 

135 y(is) = twdist(l,l,s(is)) ; 

136 end 

137 is = s > qOsolbcright; 

138 if (sum(is)) 

139 y(is) = 1; 

140 end 

141 else 

142 fprintf ( ['Don' 't know how to plot the distribution for ' 

143 ['that GSE eigenvalue. Check that the data ' ... 

144 'for it is available.']); 

145 end 

146 else 

147 fprintf ( 'Don' 't know anything about that beta!'); 

148 end 



1 global qOguess; 

2 global qOguessbclef t ; 
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3 global qOguessbcright ; 

4 global qOsol; 

5 global qOsolbcleft; 

6 global qOsolbcright ; 

7 global qlsol; 

8 global qlsolbcleft; 

9 global qlsolbcright ; 

10 load ( 'qOguess .mat ') ; 

11 qOguessbclef t = min(qOguess .x) ; 

12 qOguessbcright = max(qOguess .x) ; 

13 load CqOsol. mat') ; 

14 qOsolbcleft = min(qOsol.x) ; 

15 qOsolbcright = max(qOsol .x) ; 

16 load (' qlsol. mat') ; 

17 qlsolbcleft = min(qlsol.x) ; 

18 qlsolbcright = max (qlsol .x) ; 



1 function twmakenew 

2 global qOguessbclef t ; 

3 global qOguessbcright; 

4 global qOguess; 

5 qOguessbcright = 6; 

6 qOguessbclef t = -8; % This b.c. should never be moved 

7 "/o further to the left because of 

8 '/o singularity; function qOseed will 

9 % take care of extending range of 
10 y„ seed if need be 
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11 opts=odeset('reltol' ,le-13, 'abstol' ,le-14) ; 

12 tspan=linspace (qOguessbcright , qOguessbclef t , 1000) ; 

13 yO=[airy(qOguessbcright) ;airy(l , qOguessbcright) ;0; ... 

14 airy (qOguessbcright) "2 ; 0] ; 

15 qOguess = ode45(@deq0,tspan,y0,opts) ; 

16 save ( ' qOguess . mat ' , ' qOguess ' ) ; 

17 global qOsolbcleft; 

18 global qOsolbcright ; 

19 global qOsol; 

20 qOsolbcleft = -10; 

21 qOsolbcright = 6; 

22 options = bvpset('FJacobian' ,@qOjac, 'BCJacobian' ,@ ... 

23 qObcjac, 'reltol' ,le-13, 'abstol' ,le-14) ; 

24 grid = qOsolbcleft : (1/16) : qOsolbcright ; 

25 solinit = bvpinit(grid,@qOseed) ; 

26 qOsol = bvp4c ((9qOode,@qObc, solinit, options) ; 

27 saveCqOsol.mat' , 'qOsol') ; 

28 % $$$ Store current time/day/year in vector now 

29 now = clock; 

30 y„ $$$ Format info and put it into vector info 

31 info = [ now(4) , now(5) , ceil (now(6) ) , now(2) , now(3) , ... 

32 mod(now(l), 1000)]; 

33 % $$$ Status message prints out what stage we are in, and 

34 % $$$ what time/day/year it is. 

35 f printf ( [ 'Done with qO. Current time ' ... 

36 ['is y.02.0f:%02.0f:yo02.0f ' ... 

37 ['on day %02.0f/%02.0f/yo02.0f .\r\n'] , info); 
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38 global qlguessbclef t ; 

39 global qlguessbcright ; 

40 global qlguess; 

41 qlguessbcright = 6; 

42 qlguessbclef t = -10; 

43 opts=odeset('reltol' ,le-13, 'abstol' ,le-14) ; 

44 tspan=linspace (qlguessbcright , qlguessbclef t , 1000) ; 

45 yO= [(1/2) *airy (qlguessbcright) ; (1/2) *airy(l, ... 

46 qlguessbcright) ;0; (1/2) *airy(qlguessbcright) "2 ; 0] ; 

47 qlguess = ode45(@deql ,tspan,yO,opts) ; 

48 grid = qlguessbcleft: 1/16: qlguessbcright; 

49 global qlsolbcleft; 

50 global qlsolbcright ; 

51 global qlsol; 

52 qlsolbcleft = -10; 

53 qlsolbcright = 6; 

54 options = bvpset('FJacobian' ,@qljac, 'BCJacobian' , . . . 

55 Qqlbcjac, 'reltol' ,le-13, 'abstol' ,le-14) ; 

56 grid = qlsolbcleft : (1/16) : qlsolbcright ; 

57 solinit = bvpinit(grid,(9qlseed) ; 

58 qlsol = bvp4c (@qlode,@qlbc, solinit, options) ; 

59 save ( 'qlsol .mat 'qlsol ') ; 

60 % $$$ Store current time/day/year in vector now 

61 now = clock; 

62 '/o $$$ Format info and put it into vector info 

63 info = [ now(4) , now(5) , ceil (now(6) ) , now(2) , ... 

64 now(3), mod(now(l), 1000)]; 
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65 


y. $$$ 


Status message prints out what stage we are in, 


66 


7. $$$ 


and what time/day/year it is. 


67 


fprintf ( ['Done with ql . Current 


68 




['time is 7.02 . Of : y„02 . Of : y.02 . Of ' ... 


69 




['on day 702 . Of /%02 . Of 7X02 ' . . . 


70 




[' .Of .\r\n'] , info); 


71 


tw; 




1 


function y = matsim(beta,m,n,sigma,p, method, filename) 


2 
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7 


/. $$$ 


elements. In the case of the GSE (beta=4) n refers to 


8 


y. $$$ 


the size of the quaternion matrix; therefore the 


9 


y. $$$ 


matrices generated are then really (2n)x(2n) complex 


10 


y. $$$ 


matrices. The argument "method" refers to the method 


11 


y. $$$ 


used to generate the matrices. Setting method to "1" 


12 


y. $$$ 


uses full matrices (intuitively clearer, but also 


13 


y. $$$ 


slower) , whereas setting method to "2" uses 


14 


y. $$$ 


tridiagonal matrices a la Edelman/Dumitriu (much 


15 


y. $$$ 


faster). Also if beta is not 1, 2, or 4 then 


16 


y. $$$ 


tridiagonal matrices are used by default since they 


17 


y. $$$ 


are the only option at this time for general beta. The 


18 


y. $$$ 


eigenvalues are stored in the file "filename", the 


19 


% $$$ 


extension ".dat" is appended by default to the ouput 
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20 7o $$$ file, so an extension need not be specified. The 

21 y„ $$$ output file thus has one long column of data. The 

22 7„ $$$ first m entries are the largest eigenvalue of each of 

23 % $$$ the matrices, the next m entries are the second 

24 "/o $$$ largest eigenvalue, and so on ... 

25 7„ $$$ Record the start time 

26 begintime = clock; 

27 7o $$$ Make sure the correct number of arguments is 

28 X $$$ supplied. Else abort and exit. 

29 ifCnargin ~= 7 ) 

30 error ( ['Function should be called with seven input ' ... 

31 'arguments as '] 'matsim(beta,m,n,sigma,p, . . . 

32 method,f ilename) ' ' ; ... 

33 here ' 'beta' ' is the beta from RMT (positive . . . 

34 real), "m" . 

35 is the number of matrices in sample set (positive 

36 integer) , 

37 ' 'n' ' is the size of matrices (positive integer), 

38 '' Sigma'' is the standard deviation of normal ... 

39 distribution on diagonal matrix elements (positive 



40 



real) , 



41 



p'' is the number of eigenvalues to save ... 



42 



(positive integer), and ''fil 



ename 



is the 



43 



output file name (enclosed in single quotes; no . . 



44 



extension) . Type ' 'help matsim 



45 



to get more info.'); 



46 



elseif "isstr (filename) 
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47 error(['Last input argument must be a valid string ' ... 

48 'representing a '] 'filename']); 

49 elseif (any( [m n p method] -floor ( [m n p method])) I ... 

50 any([beta sigma]<=0)) 

51 error(['One of the arguments you gave is invalid; ' ... 

52 'type '] 'help matsim' ' to get ' 'a description ... 

53 of the arguments'' meaning and see what is wrong.']) 

54 end 

55 "/o $$$ explicitly define the square root of -1 for safety 

56 i = sqrt(-l) ; 

57 "/o $$$ Prepare output vector filled with zeros to speed up 

58 7o $$$ code by preallocating memory 

59 y = zeros (m,p) ; 

60 if (beta == 1 & method == 1) 

61 for j=l:m; 

62 A = sigma*randn(n) ; 

63 A = (A + A')/2; 

64 lambdavector = -sort (-eig(A) ) ; 

65 '/o $$$ Take the first p and store them in output vector y 

66 y(j»-) ~ (lambdavector (1 :p)) .' ; 

67 if (mod(j ,1000)==0) 

68 °/o $$$ Store current time/day/year in vector now 

69 now = clock; 

70 7o $$$ Format info and put it into vector info 

71 info = [j , now(4) , now(5) , ceil (now(6) ) , now(2) , ... 

72 now(3), mod(now(l), 1000)]; 

73 "/o $$$ Status message prints out how many matrices have been 
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74 y„ $$$ produced so far, and what time /day /year it is. 

75 fprintf ( ['\r\n Done %g runs so far. Current time ' ... 

76 ['is %g:%g:lg on ' . . . 

77 [May yog/yog/y.g.\r\n'] , info); end 

78 end 

79 elseif (beta == 2 & method == 1) 

80 for j=l:m; 

81 A = Sigma* (randn(n) + i*raiidn(n) ) ; 

82 A = (A + AO/2; 

83 lambdavector = -sort (-eig(A) ) ; 

84 y„ $$$ Take the first p and store them in output vector y 

85 y(j>-) ~ (lambdavector (1 :p)) ; 

86 if (mod(j,1000)==0) 

87 y $$$ Store current time/day/year in vector now 

88 now = clock; 

89 y $$$ Format info and put it into vector info 

90 info = [j , now(4) , now(5) , ceil (now(6) ) , now(2) , ... 

91 now(3), mod(now(l), 1000)]; 

92 y $$$ Status message prints out how many matrices have been 

93 y $$$ produced so far, and what time/day/year it is. 

94 fprintf ( ['\r\n Done yg runs so far. Current time ' ... 

95 'is y.g : yg : yg on day ' ] ' Ig/lg/tg . \r\n ' ] , inf o ) ; end 

96 end 

97 elseif (beta == 4 & method == 1) 

98 for j=l:m; 

99 X = sigma*randn(n) + i*sigma*randn(n) ; 
100 Y = sigma*randn(n) + i*sigma*randn(n) ; 
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101 A = [X Y; -conj(Y) conj(X)]; 

102 A = (A + A')/2; 

103 lambdavector = -sort (-eig(A) ) ; 

104 y„ $$$ Take the first p and store them in output vector y; 

105 % $$$ keep in min in this case eigenvalues are 

106 % $$$ doubly-degenerate, so only take every other one. 

107 y(j>:) = (lambdavector(l:2:2*p)) . ' ; 

108 if (mod(j ,1000)==0) 

109 % $$$ Store current time/day/year in vector now 

110 now = clock; 

111 % $$$ Format info and put it into vector info 

112 info = [ j , now(4) , now(5) , ceil (now(6) ) , now(2) , ... 

113 now(3), mod(now(l), 1000)]; 

114 % $$$ Status message prints out how many matrices have been 

115 % $$$ produced so far, and what time/day/year it is. 

116 fprintf ( ['\r\n Done "/og runs so far. Current time ' ... 

117 'is yog:y«g:yog on day '] ' yog/y.g/yg . \r\n , info); end 

118 end 

119 else 

120 for j=l:m; 

121 d = sigma*sqrt(chi2rnd(beta*[n:-l:l])) ' ; 

122 A = spdiags(d, l,n,n) + .. 

123 spdiags (sigma*randn(n, 1) /sqrt (2) ,0,n,n) ; 

124 A = (A + A0/sqrt(2); 

125 lambdavector = -sort (-eig(A) ) ; 

126 % $$$ Take the first p and store them in output vector y 

127 y(j.:) = (lambdavector (1 :p)) ; 
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128 if (mod(j ,1000)==0) 

129 % $$$ Store current time/day/year in vector now 

130 now = clock; 

131 % $$$ Format info and put it into vector info 

132 info = [j , now(4) , now(5) , ceil (now(6) ) , now(2) , ... 

133 now(3), mod(now(l), 1000)]; 

134 % $$$ Status message prints out how many matrices have been 

135 °/o $$$ produced so far, and what time/day/year it is. 

136 f printf ( [' \r\n Done %g runs so far. Current time ' ... 

137 ['is y.g:%g:y.g on ' . . . 

138 ['day %g/yog/y.g.\r\n'] , info); end 

139 end 

140 end 

141 y„ $$$ Done generating the eigenvalues. From this point on 

142 y„ $$$ all code is for documentation/debugging purposes. 

143 % $$$ Store current time/day/year in vector now 

144 endtime = clock; $$$ Format endtime info and put it into . . . 

145 vector info endinfo = [endtime (4), endtime (5), ... 

146 ceil (endtime (6) ) , endtime (2), ... 

147 endtime(3), endtime (1) ] ; 

148 % $$$ Format begintime info and put it into vector info 

149 begininfo = [begintime (4) , begintime(5) , ceil(begintime(6) ) , ... 

150 begintime(2) , begintime(3) , begintime ( 1) ] ; 

151 % $$$ concatenate time info into single vector 

152 timeinfo= [begininfo endinfo] ; 

153 % $$$ create a string variable that explicitly describes the 

154 y„ $$$ parameter "method" for output comments 
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155 if (method == 1) 

156 longmethod = ['full matrices (as opposed to tridiagonal ' ... 

157 'ones; see code for explanation)']; 

158 elseif (method == 2) 

159 longmethod = ['tridiagonal matrices (as opposed to full ' ... 

160 'ones; see code for explanation)']; 

161 end 

162 % $$$ Append beta as well as begin time info, n, and p to 

163 % $$$ filename 

164 filename= [filename '-' num2str (beta) '-' ... 

165 num2str(begintime(2) , '%02.0f ') ... 

166 num2str(begintime(3) , '%02.0f ') ... 

167 num2str(mod(begintime(l) ,1000) , '%02.0f ') ... 

168 num2str(begintime(4) , '%02.0f ') ... 

169 num2str(begintime(5) , '%02.0f ') '-' ... 

170 num2str(m, '7.0.0f ') num2str (n, ' -^O . Of ' ) '.dat']; 

171 % $$$ Open output file 

172 fid = f open(f ilename , ' w' ) ; 

173 % $$$ Write info blurb at the beginning of output file with 

174 % $$$ parameters. This info line is preceded by percent 

175 % $$$ sign so that it can be parsed as comment when 

176 % $$$ processing the file later to analyze the data 

177 fprintf (f id, ['%% Simulation of the first 10. Of eigenvalues ' ... 

178 ['of %0.0f matrices of size "/oO.Of by %0.0f ' ... 

179 'from the beta = "/oO.Of Hermite ensemble'] , [p ... 

180 m n n beta ] ) ; 

181 fprintf (fid, ' using /ds. ', longmethod) ; 
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182 fprintf (fid, ['Diagonal matrix elements have standard ' .. 

183 ['deviation y„f . The next /dO-Of numbers are ' 

184 ['the largest eigenvalues. The following ' 

185 ['yoO.Of numbers are the second largest eigenvalues, and ' 

186 ['so on. This simulation started at %02 . Of : %02 . Of : yo02 . Of 

187 ['on day %02 . Of /%02 . Of /yo02 . Of and ended at y„02.0f:' ... 

188 'yo02.0f :y„02.0f on day ^02 . Of /y.02 . Of /yo02 . Of . \n' ] , [sigma . 

189 m m timeinfo]); 

190 % $$$ Write eigenvalues in output file 

191 fprintf (fid, 'y.0.6f\n' ,y) ; 

192 % $$$ Close output file 

193 f close (fid); 

194 % $$$ Open temporary file for sending message 

195 fid = f open( 'matsimtempmatlabmessage . txt ' , ' w' ) ; 

196 % $$$ Put notice message in temporary file 

197 fprintf (fid, ['%% Matlab simulation is over. ' ... 

198 ['Generated the first yoO.Of eigenvalues' ... 

199 ['of yo.Of matrices of size yo.Of by y.O.Of ' 

200 ['from the beta = yo.Of Hermite ' ... 

201 ['ensemble', [p m n n beta]); 

202 fprintf (fid, 'using Xs. ', longmethod) ; 

203 fprintf (fid, ['Diagonal matrix elements have standard' . 

204 ['deviation %f . Run started at' ... 

205 ['y„02.0f :y„02.0f:y.02.0f on day' ... 

206 ['y.02.0f/y„02.0f/yo02.0f and ended at' ... 

207 [ ' yo02 . Of : y.02 . Of : y.02 .Of on day ' ... 

208 [ ' yo02 . Of /yo02 . Of /yo02 . of . ' . . . 
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209 ['\r\n' , [sigma timeinfo]); 

210 fprintf(fid, 'The data is saved in the file Zs', filename); 

211 % $$$ Close temporary file 

212 fclose(fid); 

213 % $$$ Email notice message 

214 !mail -s 'matlab message' momar@math.ucdavis.edu ... 

215 < matsimtenipmatlabmessage.txt 

216 °/o $$$ Remove temporary file 

217 !rm matsimtempmatlabmessage.txt 

218 % $$$ Exit Matlab 

219 % $$$ exit; 

1 function matsimplot (beta, p, n, m , sigma, nbins, dataf ilename , ... 

2 linespec, markersize) 

3 °/o The function matsimplot (beta, p, n, m , sigma, nbins, 

4 % dataf ilename , linespec, markersize) plots a histogram of 

5 % data generated from simulating the p~th largest eigenvalue 

6 7. of m random matrices of size nxn from a beta-ensemble 

7 7o with the program matsim. The 

8 % function expects a file named dataf ilename containing a 

9 % single column of data in which the numbers in position 

10 7o l+(p-l)*m through p*m are realisations of the p"th largest 

11 7o eigenvalue in the specified beta ensemble. The entries of 

12 7o the diagonal elements of the random matrices simulated are 

13 7o assumed to gaussian with stdev deviation sigma. In the 

14 7o beta = 1 case, the program replaces the stdev sigma on the 

15 7o diagonal by sigma/sqrt (2) which is the stdev on the off 
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16 °/„ diagonal, as it should. In the beta = 4 case it replaces 

17 % the size n of the matrices by 2n+l as it should, nbins 

18 7o number of bins. Linespec should consists of a letter 

19 % (color) and a linetype (symbol) in single quotes, with no 

20 "/o space between them. For example typing 'r-' for 

21 y„ linespecplot a red line, while 'm*' plots magenta 

22 7„ stars. Markersize gives the size (fonts) of the symbol 

23 % used to make plots 

24 if "isstr (dataf ilename) 

25 error ( ['Third input argument must be a valid string ' .. 

26 'representing an existing filename']); 

27 end 

28 data = load( ' -ascii ' , dataf ilename) ; 

29 data = data(l+(p-l)*m:p*m) ; 

30 if(beta==l I beta==2 I beta==4) 

31 if (beta == 1) 

32 Sigma = sigma/sqrt (2) ; 

33 elseif(beta == 4) 

34 n= 2*n+l; 

35 end 

36 data = ( (data*n" (1/6) )/sigma) - 2*n~(2/3); 

37 else 

38 error ( 'Don' 't know what to do for this beta!'); 

39 end 

40 [n,x] = hist (data, nbins) ; 

41 binsize = x(2)-x(l); 

42 plot(x, n/(binsize*m) , linespec , 'MarkerSize ', markersize) ; 
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hold on; 



Appendix A 



The matrix ensembles 

As mentioned in the Introduction, the three matrix ensembles discussed in this work 
were originally (and maybe more naturally) introduced as probability spaces of (self 
dual) unitary, orthogonal and symplectic matrices matrices (see ^H] for more history 
and an expanded treatment along these lines). We quickly review those set-ups here 
only insofar as they help motivate the MATLAB"""^ code we used to generate the 
eigenvalues, and set the appropriate background for the recent powerful generaliza- 
tions alluded to earlier (namely in [Hj and jH]). 

A.l The Gaussian orthogonal and unitary ensem- 
bles 

The Gaussian orthogonal (resp. unitary) ensemble GOE (resp. GUE) is defined as a 
probability space on the set Tin oi N x N real symmetric (resp. complex hermitian) 
matrices. If one requires that: 

1. the probability P{H) dH that a matrix H in 7i„ belongs to the volume element 
dH he invariant under similarity transformations H ^ AH A'' hj orthogonal 
(resp. unitary) matrices, 

127 
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2. the various Hj^ elements of the matrix H be statistically independent so that 
P{H) factors in a product of functions each depending on a single Hj^, 

then P{H) is restricted to the form 

P(i7) = Cexp(-aTrij2 + 5TriJ + c), a>0 (A.1.1) 

Converting to spectral variables and integrating out the angular part (the Vander- 
monde is the Jacobian) yields 



P^(^i,...,^^) = P^^)(^) = Cexp 



N 



ni%-^fci^ (A.1.2) 



j<k 



where /3 = 1 in the real symmetric case, and (3 = 2 in the complex hermitian case. If 
we let 6j Xj^y (3/2a + b/2a then 



Pf^ixi, ...,xn)= Pp^\x\ 



&P exp 



n . 

j<k 



Xj 



with ()A.3.3|) yielding 



(A.1.3) 



N 



(jW ^ ^2 Vr)"^/^ p-N/2-PNiN-l)/i . J-j- + 7) r(l + 2 ; 



(A.1.4) 



The following MATLAB™ code produces a GOE-distributed matrix A of size nx n 
with diagonal entries distributed A^(0,(T^), and off diagonal entries distributed as 
A^(0,(7^/2). It is adapted from [Hj; see [H], and the code listing for function matsim 
in Section 17.21 for much more efficient schemes. 

A = sigma*randn(n) ; 
A = (A + AO/2; 
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whereas the hnes below produce a GUE-distributed matrix A of size n x n also 
with diagonal entries distributed A^(0,cr^), and off diagonal entries whose real and 
imaginary parts are distributed as N{0,a^/2): 

X = sigma*randn(n) + i*sigma*randn(n) ; 
Y = sigma*randn(n) + i*sigma*randn(n) ; 
A = (A + AO/2; 



A. 2 The Gaussian symplectic ensemble 

The GSE is defined as a probability space on the set of all 2n x 2n hermitian ma- 
trices which are self-dual (in a way to be defined soon) when considered as n x n 
quaternion matrices. Arguments similar to those in the previous section lead to the 
representations ()A.1.1|) for P{H) and ()A.1.3|) for the joint probability distribution of 
eigenvalues with (3 set to 4 this time (see ^H] for details). The GSE is best defined 
in terms of quaternion matrices, since this makes its similarity to the GUE and GSE 
plain, and yields a MATLAB™ algorithm to produce the matrices. We expand on 
the construction of the ensemble here in order to motivate the MATLAB^^ code. A 
good introduction to quaternions, and the classical groups can be found in j^, which 
we use extensively in the following. Consider the algebra Q of dimension 4 over M 
generated by the basis elements Cq, ei, 62, 63 and relations 

eo Ci = 6,60 = e-i] 6^ = -60; 6i Cj = -Cj Ci = 6fc (A. 2.1) 



where 1 < i, j, k < 3 and 



I J k 



is an even permutation. Elements of Q can 



be thought of as "complex numbers with three imaginary parts" . In fact Q contains 
a sub-algebra isomorphic to the field C (the sub-algebra of quaternions of the form 
Qo + Qi Ci, each of which is mapped to the complex number qo + i qi). Therefore Q 



APPENDIX A. THE MATRIX ENSEMBLES 130 
can be considered as a vector space of dimension 2 over C. Indeed every quaternion 



g = go eo + gi ei + q2 62 + gs 63 



(A.2.2) 



can be rewritten, using the above rules, as 



g = eo(go Co + gi ei) + e2(g2 cq - gs ci] 



(A.2.3) 



which under the map to C mentioned above becomes 



g = eo(go + « gi) + e2(g2 - ^ gs) 



(A.2.4) 



Hence every quaternion can be written as a hnear combination of the two basis 
elements cq and 62 over C It is easy to check that the action of Q on itself (considered 
as a 2-dimensional C-vector space) gives a representation p : Q ^ GL{2, C) defined 
by 



p(eo) 



1 



P(ei) 



2 
-i 



(A.2.5) 



P(e2) 



v 



-1 

1 



P(e3) 



v 



-i 
-i 



This representation is often taken as the starting definition for the basis of Q. In 
this language, Q, is then thought of as an algebra (over C) of matrices with entries 
in C. We will adopt that point of view from here on (i.e. we will stop distinguishing 
between and p(ei) as we have been so far). It follows from the above discussion 
that a matrix in GL(2n,C) can be identified with one in GL{n, Q) in the following 
way; think of the complex 2n x 2n matrix as a block matrix, with blocks each 
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made up of a 2 x 2 matrices of the form 



a 



\ 



c d 



(A.2.6) 



where a, 6, c, d G C. Using the matrix representation of Q we can decompose the 
above 2x2 matrix as 



v 



a h 
c d 



J 



I i \ i 

- (a + d) eo - - (a - rf) ei + - (6 - c) 62 - - (6 + c) 63. (A.2.7) 



If each block of this form is identified with the quaternion on the right-hand side, the 
complex 2 n X 2 n matrix then becomes a quaternionic n x n matrix. We can define 
three distinct analogues of usual complex conjugation in Q. Let 



? = go eo + gi ei + q2 62 + qs 63- 



(A.2.8) 



Then following ^H] we define 



Qo eo - gi ei - q2 62 - gs 63 



go eo + g^ ei + qi. 62 + gg 63 



(A.2.9) 
(A.2.10) 
(A.2.11) 



where * denotes usual complex conjugation (recall G C). We extend this notation 
to matrices in in the usual fashion: for example if Q = {qjk} G Q" then Q = {gj^} 
and = {lij}- Let us determine what the usual matrix transformations (such as 
transposition and usual complex conjugation) of a matrix Q G GL(2n, C) correspond 
to when Q is instead considered as a matrix in Q". For example if we transpose Q 
as a complex matrix, what happens to its quaternion entry g^j? We write {Q^)ij for 
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the {i,j) quaternionic element of the matrix Q & whose transpose is taken when 
it is is considered as an element of GL{2n, C). Under transposition, the block 



a 
\^ c d J 



(A.2.12) 



in position (j, k) of Q goes to 

/ \ 

a c 

}yh d J 

in position (/c, j). The first block corresponds to the quaternion 



(A.2.13) 



1 i I i 

(Ijk = - (a + rf) eo - - (a - rf) ei + - (6 - c) 62 - - (6 + c) 63. 



(A.2.14) 



The second block is the quaternion 



kj 



1 i 1 i 

-{a + d)eo- -{a-d)ei + -{c-b)e2- -{b + c) 63 

-62 {a + d) eo + ^{a - d) ei - ^{b - c) 62 + ^{b + c) 63 ] 62 



-e2qjk 62. 



(A.2.15) 



So transposing Q as a complex matrix is the equivalent of left and right multiplying 
by 62 the quaternionic transpose of —Q. Similarly, doing this for regular hermitian 
conjugation, we see that if 



a 
d j 



\ i \ % 

- (a + rf) eo - - (a - rf) ei + - (6 - c) 62 - - (5 + c) 63, (A.2.16) 
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then 



kj 



a c 

1_- 1_- 

- (a + d) eo - - (a - d) ei + - (c - 6) 62 - - (6 + c) 63 



\ % \ % 

= - (a + rf) eo + - (a - (i) ei - - (6 - c) 62 + - (6 + c) 63 

= Q]k- (A.2.17) 
Finally we define the dual of Q G to be 

Q"" = 62 {Q') (A.2.18) 
Using the relations found above, we see that 

{Q""),^ = 62 = 62 ("62 62) 62"^ = q^,. (A.2.19) 

So this is the true analog in Q" of transpose conjugation of complex matrices. Now 
we can define the GSE properly 

Definition A. 2.1. The GSE consist of all 2n x 2n hermitian matrices which are self- 
dual when considered as n x n quaternion matrices. Equivalently, these are 2n x 2n 
hermitian matrices B that are invariant under B ^ Z B* Z^^ where Z is 2n x 2n 
block diagonal matrix with n blocks 62 down the diagonal. 

Let the elements of H be denoted by {hrs\ if H is viewed in GL{2n, C), and by 
{qjk} if it is viewed in Q" so that 1 < r, s < 2n and 1 < j,k < n. Furthermore, 
in order to be able to specify the individual quaternion components of the matrix 
elements, let us switch to the notation 

Q^k = qf^ eo + gjl^ ei + gjj + qf^ e, (A.2.20) 
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= H^ = H, 



(A.2.21; 



where we abuse notation since technically the first and second operations on H are 
in different spaces. It follows from that in terms of the quaternion entries qjk of H, 
being self-dual implies 



Qkj - Qkj - Qjk 



Combining ()A.2.22jl with ()A.2.17jl . and ()A.2.19;i yields 



(A.2.22) 



JO) _ (0) _ / (0) 
Qjk - Qkj - [Qkj 



(A.2.23) 



The last equality implies reality, and the first symmetry. Therefore if H is hermitian 
self-dual, ^ J? ^ ^ ^ forms a real symmetric sub-matrix. Similarly for 

/ = 1, 2, 3 we get 

(A.2.24) 



(0 _ _ (0 
Qjk - Qkj 



Again the last equality implying reality, and the first antisymmetry. Therefore if H 
is hermitian self-dual, l^j^ | ^ j^k < n forms a real antisymmetric sub- matrix for 
each /. Finally we need to know how to go back from the qj2 to the hrs- It is easy to 
see that if 



a 
y c d J 



(A.2.25) 



corresponds to 



1 i 1 i 

Qjk = ^{a + d)eo- -{a-d)ei + -{b-c)e2- -(b + c) 63, 



(A.2.26) 
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then 



(0) , (1) , (2) , (3) 



(A.2.27) 



corresponds to 



with 



a 
\^ c d J 



(A.2.28) 



a 
b 
c 
d 



(0) , • (1) 

q-jk + ^ q)k 

(2) , ■ (3) 

qjk + « q)k 

, (2) • (3) 

- 1 q)k - ^ q)k 

(0) ■ (1) 

- ^ qjk 



(A.2.29) 
(A.2.30) 
(A.2.31) 
(A.2.32) 
(A.2.33) 



which yields the following recipe to construct a random 2n x 2n GSE matrix M: 

1. generate a random n x n real symmetric matrix and 3 random n x n real 
antisymmetric matrices X, Y, Z. 

2. the sub- matrix A of M which in our notation corresponds to the entry a in each 
complex 2x2 block will be given hj A = W + iX. Similarly B = Y + i Z, 

C = -{Y -iZ), D = W -iX 

3. then M = A(g)Ii + B(g)l2 + C(S)l3 + D(g)h where 



Ol,k 02,k 
<^3,fc <^4,fc 



(A.2.34) 



This recipe is simplified by noting that the matrix M above has the same eigenvalues 
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as the matrix 



H 



\ 



C D 



(A.2.35) 



Note that if A denotes the complex conjugate of A, and A* its transpose, then by 
construction A = A^, A = D, B* = —B, —B = C . Hence we can equivalently let 



H 



( A B^ 



\ 



-B A 



(A.2.36) 



where A and B are complex matrices satisfy = A and 5* = —5. Hence the 
following MATLAB™ code produces a GSE-distributed matrix A of size 2n^1n 
with diagonal entries distributed iV(0,(T^), and off diagonal entries distributed as 
A^(0,(T^/2). It is adapted from [Hj; see [H], and the code listing for function matsim 
in Section 17.21 for much more efficient schemes. 

X = sigma*randn(n) + i*sigma*randn(n) ; 
Y = sigma*randn(n) + i*sigma*randn(n) ; 
A = [X Y; -conj(Y) conj(X)]; 
A = (A + AO/2; 

A. 3 Selberg's integral 

For any positive integer n\%i dx = dx\ - ■ ■ dxn and define 

A(x) = A(xi, . . . ,x„) = JJ {xj-Xk), 

l<j<k<n 

if n > 1, and A(a;) = 1 if n = 1. Also let 



$(x) = ...,x„)= I A(x)i'^ n s""' (1 - ^^•)''"'- 
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Then, 

Theorem A. 3.1 (Selberg, Aomoto). 



„1 „i n-l 

I{a,p,-f,n)= ■■•/ <!){x)dx = T\ 

Jo Jo „_n 



r(l + 7 + j7)r(a+j7)r(/3 + j7) 
^.^^ r(l+7)r(a + /3+(n + j-l)7) ' 



(A.3.1) 

and for 1 < m < n, 

■■■ f x,X2---Xm ^x)dx = TT "^/j" ~ ^ n / • • • / 

(A.3.2) 

valid for integer n, and complex a, (3, 7 with Rea > 0, Re/? > 0, and 



Re 7 > — min 



1 Re a Re /? 
n' n — 1 ' n — 1 



The identity ()A.3.1|) is the one commonly referred to as Selberg's integral, and 
()A.3.2|) is Aomoto's extension of it. See ^H] for proofs and an expanded discussion 
of its consequences and applications. The following specialization is of particular 
interest here; put Xj = yj/L, a = P = a L"^ + 1, and take the limit L — > 00 in ()A.3.1jl 
to obtain 



00 poo -- 

|A(x)|^^ ]^exp(— ax^) dx 



00 J -00 



= (2 (2 . n r(i+Vr{^+,) - 

Set 7 = /3/2 and a = 1 to justify (fTX^ . whereas a = {3/2 gives (IA.1.4jl . 
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As mentioned in Section ESI the choice of standard deviation is not the same in all 
three ensembles. This affects the way we scale data obtained in MATLAB™ simu- 
lations with the code give in the previous sections in this appendix. We treated the 
/? = 2,4 cases with weight function exp(— x^) corresponding to a standard deviation 
(Td = 1/ \/2 on the diagonal matrix elements, and the scaling 



V2iV- 



Therefore if we want to rescale data with an arbitrary standard deviation on 
diagonal elements, we need to let t — > tj^^a^ resulting in the more general change 
of variables 

In the /5 = 1 case however our choice of weight function exp(— a;^/2) corresponds to 
(Td = 1. A similar argument would thus yield 

In order to preserve some uniformity in the formulas, we can use the standard devia- 
tion of the off diagonal elements instead in the /3 = 1 case. Indeed, if we denote by cTq 
the standard deviation on the off-diagonal elements then Od = by construction 

in all three ensembles. Hence in our treatment, an = I/V2 in the /? = 1 case, and for 
general (Tq we have, similar to the /5 = 2,4 change of variables 
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